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Abstract 

In  the  framework  of  the  proposed  "continuous  approach"  to  constrai¬ 
ned  optimization  problems,  we  describe  two  new  solution  methods  which  re¬ 
sulted  from  the  research.  The  first  is  a  continuous  "inexact"  method  for  sol¬ 
ving  systems  of  nonlinear  equations  and  complementarity  problems  (along 
the  lines  of  the  DAFNE  Method),  and  the  second  is  a  continuous  method 
for  solving  the  linear  programming  problems  (along  the  lines  of  Karmarkar*s 
method)  which  is  shown  to  be  quadratically  convergent. 

Some  numerical  experience  on  a  number  of  test  problems  is  reported. 
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1.  Introduction 

This  is  the  final  report  on  the  work  performed  from  September  1986 
to  December  1991,  under  contract  n.  DAJA  45-86-C-0028  awarded  to  the 
University  of  Rome  "La  Sapienza”  on  the  research  project  "Numerical  Opti¬ 
mization",  by  the  principal  investigator  Francesco  Zirilli  and  his  co-workers. 

The  objective  of  the  research  is  described  in  par.  2,  the  results  of  the 
research  are  described  in  par.  3,  and  some  conclusions  are  in  par.  4. 


2.  Objective  of  the  research 

The  subject  of  the  research  was  the  field  of  those  problems  in  con¬ 
strained  optimization  which,  starting  from  the  linear  programming  problem, 
can  be  formulated,  with  growing  degree  of  generalization,  first  as  linear 
complementarity  problems  and  second  as  nonlinear  complementarity  pro¬ 
blems. 

The  objective  of  the  research  was  to  attack  the  above  problems  by 
means  of  the  so-called  "continuous  approach"  to  optimization  (as  opposed  to 
the  so-called  "pivotal"  methods,  such  as  the  simplex  method  for  linear  pro¬ 
gramming),  with  special  consideration  for  the  interesting  cases  of  non-con- 
vex  or  ill-conditioned  problems,  and  problems  with  a  very  large  number  of 
variables;  and  in  particular  the  objective  was  to  investigate  the  possibility  of 
applying  to  the  above  problems,  suitably  transformed  into  nonlinear  equa¬ 
tion  problems,  the  methods  developed  by  the  prindpal  investigator  and  his 
co-workers  for  solving  nonlinear  equations  and  global  optimization  pro¬ 
blems,  based  on  the  numerical  integration  of  suitable  ordinary  or  stochastic 
differential  equations  (refs.  [1]  to  (4]). 


3.  Results  of  the  research 


During  the  development  of  the  research  the  complementarity  pro- 
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blems  proved  to  be  much  more  difficult  than  it  had  been  anticipated  and  it 
became  clear  that  the  original  plans  where  somehow  too  ambitious. 

The  final  outcome  of  the  research,  if  judged  against  the  original 
plans,  is  therefore  admittedly  less  satisfactory  than  it  was  originally  hoped; 
nevertheless  a  number  of  interesting  results  have  been  obtained,  so  that  we 
feel  that  the  research  is  still  to  be  considered  at  least  partially  successful. 

The  main  results  of  the  research  are  contained  in  the  two  papers 
numbered  [5]  and  [6]  in  the  list  of  references,  which  are  described  in  the  fol¬ 
lowing  paragraphs  3.1  and  32,  and  enclosed  as  Appendix  1  and  Appendix  2. 

Report  on  some  work  performed  in  other  directions  along  the  lines  of 
the  original  research  plan,  together  with  some  related  results  of  auxiliary 
and  preliminary  nature,  were  described  in  the  Periodic  Technical  Reports; 
see  also  the  papers  numbered  [8]  and  [9]  in  the  list  of  references. 

The  research  has  also  stimulated  scientific  contacts  with  several  Ita¬ 
lian  and  foreign  scholars. 

The  above  results  have  been  disseminated  by  means  of  the  aforemen¬ 
tioned  papers  on  high-standard  academic  journals,  and  seminars  at  Accade- 
mia  dei  Lincei,  Rome  (ref.  [7],  which  originated  paper  [8]),  at  two  meetings 
of  CECAM,  Centre  Europ^en  de  Calcul  Atomique  et  Mol^culaire,  the  first 
in  Ermelo  (The  Netherlands),  ref.  [10],  and  the  second  at  CECAM  main  of¬ 
fice  in  Orsay  (Paris),  France,  ref.  [11]. 


3.1.  The  first  paper 

The  first  paper  (ref.  [5],  and  Appendix  1)  can  be  summarized  as  fol¬ 
lows. 

A  class  of  algorithms  is  developed  for  the  numerical  solution  of  non¬ 
linear  systems  of  equations  and  complementarity  problems,  based  on  the 
fact  that  the  solution  of  complementarity  problems  can  be  reduced  to  the  so¬ 
lution  of  systems  of  nonlinear  equations  by  means  of  a  transformation  first 
suggested  by  Mangasarian. 

The  method  is  "continuous"  since  it  looks  for  the  solution  of  the  non- 
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linear  system  by  following  the  numerical  solution  trajectories  of  a  suitable 
differential  equation,  as  in  previous  work  of  the  same  authors  such  as  the 
method  implemented  in  the  package  DAFNE,  described  in  Refs.  [1]  and  [2]. 

At  each  numerical  integration  step,  the  DAFNE  method  requires  the 
solution  of  an  NxN  system  of  linear  equations,  and  the  cost  of  solving  such  a 
system  when  a  large  numer  N  of  unknowns  is  involved  is  the  most  important 
part  of  the  computation. 

The  present  method  can  be  called  "inexact",  since  it  computes  only  an 
approximate  solution  of  the  above  linear  system,  by  means  of  a  conjugate- 
gradient  procedure  which  is  suitably  stopped  before  "convergence",  i.e.  after 
a  number  mC<N)  of  steps  depending  on  the  norm  of  the  residual.  For  these 
algorithms  local  convergence  and  Q-superlinear  rate  of  convergence  has 
been  proved.  The  algorithms  have  been  used  to  solve  three  complementarity 
problems  derived  from  variational  inequalities  of  mathematical  physics  very 
successfully.  The  complementarity  problems  considered  had  up  to  900  varia¬ 
bles. 


32.  The  second  paper 

The  second  paper  (ref.  [6]  and  Appendix  2)  can  be  summarized  as 
follows. 

The  paper  introduces  a  new  method  for  solving  the  linear  program¬ 
ming  problem,  i.e.  the  problem  of  minimizing  a  linear  cost  function  of  seve¬ 
ral  real  variables,  subject  to  linear  equality  and  inequality  constraints. 

Following  Karmarkar  [12],  the  paper  considers  the  problem  in  the 
"canonical"  form 

minimize  f(x)  =  f^x 

subject  to 

Ax 

iTx  =  1 

xj^O,  i=l, ....  n 


-  6 


where 


=  (1.U.1)’ 


are  real  column  vectors  with  n  elements,  A  is  a  real  m  x  n  matrix  of  rank  m, 
with  Ag=i),  n^2,m<n,  and  without  loss  of  generality  the  objective  fun¬ 
ction  f(x)  may  be  "normalized",  i.e.  fOf*)=0  if  is  a  solution  of  the  problem. 

In  this  pq)er  it  is  shown  that  Karmarkar’s  method  [12]  is  in  fact  equi¬ 
valent  to  applying,  to  a  suitable  initial  value  problem  for  a  system  of  ordina¬ 
ry  differential  equations,  the  numerical  integration  method  known  as  Euler’s 
method  with  variable  stepsize,  and  obtaining  the  problem  solution  ic*  as  the 
limit,  as  t  goes  to  infinity,  of  the  numerically  computed  solution  x(t)  to  the 
initial  value  problem,  starting  from  the  initial  point  =(l/n)s. 


The  proposed  method  is  also  based  on  the  above  interpretation  of 
Karmarkar’s  method,  but  with  two  main  differences: 


1)  the  initial  value  problem  is  based  on  a  different  system  of  ordinary  diffe¬ 
rential  equations; 

2)  the  numerical  integration  method  is  a  linearly  implicit  A-stable  method 
with  variable  stepsize. 

The  resulting  algorithm  is  shown  to  be  quadratically  convergent. 

The  computational  cost  of  one  step  of  the  proposed  algorithm  is 
shown  to  be  of  the  same  order  of  one  step  of  Karmarkar’s  algorithm. 

While  one  step  of  the  classical  simplex  algorithm  [13]  for  linear  pro¬ 
gramming  is  much  cheaper,  it  may  be  expected  that  -  due  to  the  quadratic 
convergence  -  the  number  of  iterations  needed  to  solve  a  linear  program¬ 
ming  problem  to  a  given  accuracy,  should  be  approximately  independent  of 
the  problem  size  n. 

Some  numerical  results  are  also  reported,  which  appear  to  support 
such  expectatioa 

The  algorithm  was  tested  on  ten  test  problems,  one  originating  from 
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the  operations  of  an  industrial  plant  in  central  Italy,  and  the  other  nine  pro¬ 
vided  by  the  System  Optimization  Laboratory  at  Stanford  University. 

The  results  are  reported  in  Table  1  where  n  is  the  number  of  varia¬ 
bles,  m  the  number  of  constraints,  k  is  the  index  of  the  first  step  that  verifies 
the  stopping  rule 

f(j)<  10-*.f(jo) 


and  is  the  corresponding  value  of  f(x). 

We  note  that  the  test  problems  with  n,m  ^  5  are  solved  in  about  ten 
steps,  and  that,  while  n  and  m  vary  by  an  order  of  magnitude,  the  number  k 
of  steps  needed  to  solve  the  problem  varies  only  by  a  factor  of  two. 

TABLE  1 


Test  problem 

m 

n 

k 

‘'k 

1.  ZIRl 

304 

543 

21 

2.41D-10 

2.  ADLTTTLE 

57 

141 

21 

3.16D-09 

3.  AHRO 

28 

54 

12 

1.52D-12 

4.  BEACX)NFD 

173 

298 

20 

3.91D-09 

5.  BLEND 

75 

117 

21 

1.47D-12 

6.  ISRAEL 

175 

319 

17 

1.94D-10 

7.  SC105 

106 

166 

13 

136D-11 

8.  SC50A 

51 

81 

14 

1.48D-14 

9.  SC50B 

51 

81 

11 

7.84D-10 

10.  SHARE2B 

97 

167 

21 

1.78D-10 

4.  Conclasions 


The  research  resulted  in  two  new  methods,  one  for  solving  comple¬ 
mentarity  problems,  and  the  other  for  solving  the  linear  programming  pro¬ 
blems,  both  based  on  the  so-called  "continuous"  approach  to  optimization. 
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Successful  solution  was  obtained  for  the  complementarity  problems 
on  test  problems  with  up  to  900  variables. 

However,  due  to  the  great  difficulty  of  complementarity  problems  -  in 
fact  much  greater  than  expected  -  the  hoped  for  attack  of  much  more  diffi¬ 
cult  problems  proved  to  be  unsuccessful. 

The  linear  programming  method  was  shown  to  be  quadratically  con¬ 
vergent,  and  was  successfully  tested  on  preliminary  test  problems  with  up  to 
about  300  variables  and  540  constraints. 
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An  Inexact  Continuons  Method  for  the  Solution 
of  Large  Systems  of  Equations  and 
Complementarity  Problems 

F.  ALUFFI-PENTINI  •  V.  PAJUSI  -  F.  ZIRILLI^*) 


Dedkito  alia  memofia  di  Carlo  Cattaneo.  maestro  ed  amko 


Rl  ASS  UNTO  -  Si  considera  un  nuovo  metodo  per  la  risoliuione  rtumerica  tia  di 
sistemi  di  eqvazioni  nan  lineari  sia  di  ptvblemi  di  eomplementaritd,  che  si  basa  sul  fatto 
che  la  risoluzione  di  problemi  di  eomplementaritd  si  pud  ricondurre  alia  risoluzione  di 
sistemi  di  ejuazioni  non  lineari  mediante  una  trasformazione  suggerila  da  ifangasar- 
ion.  n  metodo  e  'eontinuo"  in  quanto  la  soluzione  del  sistema  viene  cercata  seguendo 
le  traieltorie  ottenute  per  integrazione  numerioa  di  un'opportuna  equazione  differenziale 
—  come  in  precedenti  lavori  degli  autori  —  e  st  pud  dire  ‘'inesalto"  nel  senso  che  fa 
uso  di  un  metodo  di  gradienti  coniugati,  opportunamente  arrestato  "prima  della  con- 
vergerua",  per  la  risoluzione  del  sistema  lineare  che  nasoe  nell 'integrazione  numerica 
dell 'equazione  differenziale.  B  metodo  appare  partieolarmente  efficiente  per  problemi  in 
cui  compare  un  gran  numero  di  variabili  indipenderUi,  nei  quali  la  parte  prevalente  dello 
sforzo  di  calcolo  e  rappresentata  dalla  soluzione  di  un  sistema  lineare  ad  ogni  pasao  di 
integrazione.  Vengona  dimostrate  la  convergenza  locale  e  la  convergenza  Q-supertineare 
del  metodo,  e  vengona  presentati  olcuni  risultati  numeric!  relativi  o  problemi  di  com~ 
plementaritd  della  fisioa  matematioa. 

Abstract  -  H^e  coneider  a  new  method  for  the  numerical  solution  both  of  non¬ 
linear  systems  of  equations  and  of  complementarity  problems,  based  on  the  fact  that 


^*^The  research  reported  in  this  document  has  been  made  possible  through  the  support 
and  sponsorship  of  the  U.S.  Government  through  its  European  Research  OfRce  of  the 
U.S.  Army  under  contract  n.  DAJA  45-86-C-0028. 
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F.  ALUFFI-PENTINI  -  V.  PARISI  -  F.  ZUULU 


[2] 


the  solution  of  complementarity  problems  can  be  reduced  to  the  solution  of  nonlinear 
systems  of  equations  by  means  of  a  transformation  first  suggested  by  Afangasanan.  The 
method  is  ‘continuous’  since  it  looks  for  a  solution  of  the  nonlinear  system  by  following 
the  numerical  solution  trajectories  of  a  suitable  differential  equation  —  as  in  previous 
work  by  the  present  audtors  —  and  can  be  called  ‘inexact’  since  it  uses  a  conjugate- 
gradient  method  which  is  suitably  stopped  "before  convergence’  for  the  solution  of  the 
linear  system  arising  in  the  numerical  integration  of  the  differential  equation.  The 
method  appears  to  be  particularly  effective  for  problems  involving  a  large  number  of 
independent  variables,  where  the  computational  cost  is  dominated  by  the  solution  of  a 
linear  system  at  each  integration  step.  Local  convergence  andQ  :-j:~’-':near  convergence 
of  the  method  are  proved,  under  suitable  assumptions,  and  some  numencal  experience 
on  complementarity  problems  of  mathematical  physics  is  presented. 

Key  Words  -  Numerical  analysis  -  Nonlinear  equations  -  Mathematical  pro¬ 
gramming  -  Complementarity  problems. 

A.M.S.  Classification:  65H10  -  65K05 


1  -  Introduction 

Let  be  the  TV^-dimensional  real  enclidean  space,  let  x  = 
{xi,X2,..,.,xnY  €  IR^  be  a  vector,  and  for  x,y  6  IR‘'^  let  (x,y)  = 

N 

l|x||  =  be  the  euclidean  scalar  product  and  norm;  where 

1=1 

necessary  ||  •  ||  will  indicate  also  the  matrix  norm  induced  by  the  eu¬ 
clidean  vector  norm.  Given  f  :  IR^  —*■  IR^  we  will  be  concerned  with  two 
classes  of  problems  in  this  paper:  the  problem  of  solving  the  system  of 
simultaneous  nonlinear  equations 

(1.1)  f(x)  =  0 

that  is:  find  x*  6  IR^  such  that  f(x’)  =  0,  and  the  complementarity 
problem 

(1.2)  X  >  0 

(1.3)  f(x)  >  0 

(1.4)  (x,f(x))  =  0 

where  x  >  0  means  Xj  >  0,  i  =  1,2, . . . ,  and  similarly  f(x)  >  0  means 
/•(x)  >  0,  i  =  1,2,..., IV,  /i(x)  being  the  components  of  f,  that  is:  find 
X*  such  that:  x*  >  0,  f(x’)  >  0,  (x*,f(x*))  =  0. 


(3] 


An  Inexact  Continooos  Method  for  the  Solation  etc. 
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The  importance  of  the  problem  of  solving  a  system  of  simultaneous 
equations  is  well  known.  When  f(x)  =  Ax  +  b  is  an  affine  map  the 
(linear)  complementarity  problem  has  been  considered  by  Cottle  and 
Dantzig  in  [1]  and  contains  as  special  cases  the  linear  prograunming  and 
the  quadratic  programming  problem.  In  the  case  when  f(x)  is  a  possibly 
nonlinear  function  of  x  the  (nonlinear)  complementarity  problem  is  a 
rather  general  problem  and  contains  as  special  cases  the  Kuhn- Tucker 
first-order  necessary  conditions  for  the  nonlinear  programming  problem 
and  has  been  widely  studied,-  see  for  example  GOULO  and  TOLLE  [2]. 

The  linear  and  nonlinear  complementarity  problems  have  applica¬ 
tions  in  such  diverse  areas  of  flow  in  porous  media  [3],  image  reconstruc¬ 
tion  [4],  [5],  game  theory  [6]- 

In  this  paper  we  will  be  concerned  with  the  problem  of  the  numerical 
solution  of  nonlinear  systems  of  equations  and  complementarity  problems. 
Usually  complementarity  problems  are  approached  numerically  with  piv¬ 
otal  methods  (for  example  the  simplex  method  for  linear  programming). 
The  pivotal  methods  are  usually  of  the  “step  by  step”  improvement  type, 
that  is,  given  a  problem  for  which  a  solution  is  sought,  the  standard 
approach  is  to  attempt  to  define  recursively  a  sequence  of  approximate 
solutions  which  have  the  basic  property  of  making  an  improvement  in  a 
suitable  “objective  function”.  When  the  problem  satisfies  some  convexity 
and/or  monotonicity  assumptions  the  pivotal  methods  are  guaranteed  to 
converge  and  if  only  a  moderate  number  of  independent  variable  is  in¬ 
volved  (up  to  few  hundreds)  their  numerical  performance  is  satisfactory. 

In  recent  years  there  has  been  a  growing  interest  in  the  use  of  con¬ 
tinuous  methods  in  nonlinear  optimization;  see  for  example  Allgower 
and  Georg  [7]  for  a  review  of  simplicial  methods  in  the  computation 
of  fixed  points  and  the  solution  of  nonlinear  equations,  and  Bayer  and 
Lagarias  [8]  for  the  interpretation  of  Karmarkar’s  linear  programming 
algorithm  as  a  method  that  follows  a  trajectory  of  a  suitable  system  of 
ordinary  differential  equations.  In  particular  the  present  authors  have 
developed  a  method  for  solving  systems  of  nonlinear  equations  bared  on 
the  numerical  integration  of  an  initial-value  problem  for  a  system  of  or¬ 
dinary  differential  equations  inspired  by  classical  mechanics  [9],  [10],  [11], 
[12]  and  a  method  for  global  optimization  based  on  the  numerical  inte¬ 
gration  of  an  initial  value  problem  for  a  system  of  stochastic  differential 
equations  inspired  by  statistical  mechanics  [13],  [14],  [15].  In  section  2  the 
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algorithms  introduced  in  [10]  to  solve  systems  of  nonlinear  equations  are 
modified  to  obtain  an  “inexact”  solution  of  the  linear  systems  appearing 
in  each  iteration  in  the  spirit  of  Dembo,  Eisenstat  and  Steihaug  [16]. 
These  new  algorithms  are  particularly  effective  for  problems  involving  a 
large  number  of  independent  variables  where  the  computational  cost  is 
dominated  by  the  solution  of  the  linear  system  at  each  step.  Under  suit¬ 
able  hypotheses  local  convergence  and  Q-superlinear  convergence  of  these 
new  “inexact”  algorithm  for  nonlinear  systems  of  equations  are  proved. 
In  section  3  the  complementarity  problem  is  transformed  into  a  nonlin¬ 
ear  system  of  equations  following  Mangasarian  [17]  and  therefore  the 
algorithms  previously  developed  provide  a  class  of  locally  convergent  Q- 
superlinear  methods,  which  are  not  of  the  “step-by-step  impro\ement” 
type,  for  the  solution  of  complementarity  problems.  Finally  in  section  4 
some  numerical  experience  obtained  with  the  algorithms  of  section  2  and 
3  on  some  complementarity  problems  of  mathematical  physics  is  shown. 

Some  of  the  results  of  this  paper  have  been  announced  in  [18]. 


2-  Some  inexact  algorithms  for  nonlinear  systems  of  equations 

Let  f(x)=(/i(x),/7(x), . . .  ,/Af(x))’’  e  where  /^(x),  i=  1,2, . . . , 
are  real-valued  regular  fimctions  defined  for  x  =  (ii,  zj, . . . ,  Xff)^  G  IR^. 
In  order  to  solve  the  system  of  simultaneous  equations 

(2.1)  f(x)  =  0 
we  define 

(2.2)  /•(*)= 

i=l 

It  is  easy  to  see  that  x*  is  an  isolated  minimizer  of  F{x)  and  F(x*)  = 

0. 

In  [9],  [10],  [11],  [12]  the  idea  has  been  proposed  and  developed  of 
associating  to  the  nonlinear  system  (2.1)  the  following  system  of  second- 
order  ordinary  differential  equations: 

A'T 

(2.3)  »‘|^(‘)  =  f  (<)  -  Vf  (x(l))  «  e  10,  +00) 
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Where  D  is  a  AT  x  JV  positive  symmetric  matrix,  n,g  are  positive  con¬ 
stants,  VF(x)  is  the  gradient  of  the  function  F(x)  with  respect  to  x. 
The  equation  (2.3)  represents  Newton’s  second  law  (mass  x  acceleration 
=  force)  for  a  particle  of  mass  moving  in  IR^  subject  to  the  force  -VF 
given  by  the  potential  F  and  to  the  dissipative  force  -gDdxfdt. 

If  X*  is  an  isolated  minimizer  of  F(x)  then  x(t)  =  x*,  V  t  6  [0,-l-oo), 
is  a  solution  of  (2.3);  consider  the  Cauchy  data 

(2.4)  x(0)  =  ^o 

(2.5)  ^(0)  =  1^0 

and  let  x(t,^o,fJo)  be  the  solution  of  the  Cauchy  problem  (2.3),  (2.4), 

(2.5) . 

It  can  be  shown  that  there  exists  a  neighborhood  U  C  of  |**|  6 


IR*^  such  that  if 

®  6  £/’  we  have: 

0 

(2.6) 

Urn  ||x(t,^o,*lo)-x*|l  =  0 

C—»00 

Hence  in  order  to  solve  the  system  of  nonlinear  simultaneous  equa¬ 
tions  by  integrating  numerically  the  Cauchy  problem  (2.3),  (2.4),  (2.5), 
we  are  primarily  interested  in  the  equilibrium  points  reached  asymptoti¬ 
cally  by  the  trajectories  of  (2.3)  (hopefully  solutions  of  (2.1))  rather  than 
in  the  accuracy  of  the  numerical  scheme.  So  that  of  particular  interest  are 
numerical  methods  enjoying  a  special  stability  property  called  i4-stability 
(lO). 

Let  /  6  IR,  let  y,(g  €  IR"*  and  ^(t,y)  6  IR”*  be  a  given  function 
continuous  in  t  and  continuously  differentiable  with  respect  to  y,  such 
that  the  initial-value  problem: 

(2.7)  ^(0  =  vK«.y)  f€(0,-|-oo) 

(2.8)  y(0)  =  ^0 


has  a  solution  y(t,fo)  for  t  6  [0,-boo). 
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The  simplest  choice  of  i4-stable  linearly  implicit  method  to  integrate 
numerically  (2.7),  (2.8)  is: 


(2.9)  (/-ft*n)(y«+i -y«)  =  hv>n  n  =  0,1,2,... 

(2.10)  yo  =  fo 


where  y„  is  the  numerically  computed  approximation  of  y(n/i,fo))  ^  is 
the  identity  matrix  acting  on  IR”,  h  >  0  is  the  stepsize,  for  n  =  0,1,2,. . . 
i„  =  n/i,  =  v(t„,y„),  in  =  iitn,yn)  where  $(t,y)  =  dip/dy  is  the 
Jacobian  of  ip  with  respect  to  y.  We  note  that  when  ip{t,y)  =  Ay  is  a 
linear  map  (2.9)  reduces  to  the  backward  Euler  method. 

.After  rewriting  (2.3)  as  a  first-order  system 


(2.11) 

(2.12) 


dt 

dv 

It 


-  I Dyr-  -VF{x) 

p  p 


formulae  (2.9),  (2.10)  with  variable  stepsize  A„,n  =  0,1,...  (i.e.  to  =  0, 
t„  =  X)  ^*1  ”  =  li2, ...)  are  applied  to  (2.11),  (2.12),  (2.4),  (2.5).  In  this 

i=0 

case  the  map  <p :  IR*^  — ►  IR*^  will  be  given  by 


(2.13) 


V 

-i£»v-iV/’(x) 


so  that  its  Jacobian  matrix  is  given  by 


(2.14) 


where 


I 


(2.15) 


I(x)  =  2 


J{xfJ{x)  +  Y,fi(x)Hi{x) 

>=i 


J{x)  =  d{{x)f&x  is  the  Jacobian  of  f  with  respect  to  x  and  n,{x)  is  the 
hessian  of  /i(x). 
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[]] 

Let  s„  =  x„+i  -  x„,  n  =  0,1,2,...;  after  some  simple  algebra  (2.9) 
becomes: 


(2.17) 

v„+i  =  ^  n  =  0,1,2,... 

(2.18) 

x„+l  =X„+8, 

where  I„  =  I(x„),  Vfn  =  VF(Xn).  In  order  to  avoid  the  computation 

of  Hi{x),  i  =  1,2, . . . ,  iV,  at  each  iteration  and  since  we  are  looking  for 

s 

points  x*  such  that  f(x*)  =  0  the  term  fi{x)  in  (2.15)  is  dropped  so 

•  si 

that  L{x)  is  substituted  by 

(2.20) 

i(x)  =  2J^(x)J(x). 

Equation  (2.16)  will  be  replaced  by 

(2.21) 

where 

(2.22) 

A(x,h)  =  i(x)  +  i[^/+«7Z)] 

and 

(2.23) 

An  —  ^(Xn,  h„) 

(2.24) 

b„  =  -VFn  + 

we  note  that  the  matrix  An  is  symmetric  and  positive  definite. 
We  have  the  following  theorem: 
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Theorem  2.1.  Let  f  :  be  tvnce  continuotisly  differen¬ 

tiable,  F(x)  =  f(x)^f(x)  and  L{x)  be  given  by  (2.15).  Let  x*  6  be 
such  that  f(x*)  =  0,  J(x’)  is  nonsingular  (i.e.  x‘  is  a  nondegenerate 
solution  of  the  system  (2.1)j  and  the  following  Lipschitz  conditions  holds: 

(2.25)  llX(x)  -  Z(x-)||  <  7||x  -  x*|l  V  x  6  5  =  {x|  l|x  -  x*)!  < 

for  some  constants  7  and  S  greater  than  zero.  In  the  iteration  (2.21), 
(2.17),  (2.18)  let  {/»„},  n  =  0,1,2,...,  be  a  sequence  of  positive  numbers 
such  that 

(2.26)  lim  hn  =  00 

n— *00 

then  there  exists  h  >  Q  such  that  for  hn  >  h,  n  =  0, 1, . . . ,  x*  is  a  point 
of  attraction  o/(2.21),  (2.17),  (2.18)  and  the  rate  of  convergence  is 

(i)  Q-superlinear  if  h';'’  <  7i|lVF(x„)l|,  71  >  0,  n  >  tiq,  for  some  71, 

no  >  0. 

(iij  Q-quadratic  if  h-^  <  72||Vf’(x„)ll’,  7j  >  0,  n  >  no,  for  some  7j, 

no  >  0. 

Proof.  Let  us  rewrite  (2.21),  (2.17),  (2.18)  as 

(2.27)  x„+i  =  G(x„,/i„)+ r— — A;‘(x„-Xn_i)  n  =  0,1,2,... 

nn^n-l 


where 

(2.28)  G(x,/i)  =  X  -  >l(x,h)-‘V/’(x) 

with  the  initial  conditions  xo  =  (o>X-i  =  (0  ~  ^-1  = 

that  is  (2.21),  (2.17),  (2.18)  can  be  interpreted  as  a  two-step  iteration. 
Since  X*  is  a  nondegenerate  solution  of  the  system  (2.1)  x*  is  an  isolated 
minimizer  of  F(x)  and  Vf’(x*)  =  0.  Moreover  for  h  >  0  the  symmetric 
matrix  A{x,h)  is  positive  definite  so  that  A(x,/t)~‘  exists  that  is  G{x,h) 
is  well  defined  for  x  €  IR"'^  and  h  >  0  and  x’  is  fixed  point  of  G{x,h). 
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Let  /3  =  l|L(x*)~'|l  and  let  e  6  (0,(2/3)"^)  then  there  exists  S  >  0 
and  h  >  0  such  that: 


(2.29) 


l|I(x’)  -  A(x,/i)||  <  £  Vx  6  S  =  {x|  llx  -  x*ll  <  «} 

Vh  >  /» 


In  fact 


l|I(x-)  -  A{x,h)\\  <  l|L(x-)  -  i(x)ll  +  lll(x)  -  A(x,/»)11 
since  I(x*)  =  i(x*)  there  exists  6  such  that: 

\\L{x')-i{x)\\<^£  Vx6S 
and  for  a  suitable  h>  0 

(2.30)  lll(x)  -  A{x,h)\\  =  ^11^^  +  9D\\  <\s  V  h  >  h 

From  (2.29)  and  the  perturbation  lemma  (lemma  2.3.2  p.45  of  Or¬ 
tega  and  Rheinboldt  [19])  it  follows  that  A(x,/i)“’‘  satisfies 

(2.31)  ||A(x,/i)-‘l|  <  a  =  VxeS,V/i>h 

Moreover 


(2.32)  ||G(x,h)-x-||  <u(x./t)||x-x*||  VxeS,V/»>h 
where 

(2.33)  u(x,h)  =  a[||A(x,h)  -  X(x)ll  P(x)  -  i(x*)l|  -I-  l|q(x)||] 
and 


fO 


?(x)  =  { 


[|Vf(x)- VF(x-)  -  L{xn(x  -  x*)|| 


X-X* 


X=  X* 

xjix' 
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In  fact 

\\G{x,h)  -  x’ll  =  ||A(x,/»)-‘  [a(x,/i)(x  -  X*)  -  Vf  (x)]  |1  < 

<  a{[||A(x,/i)  -  I(x)||  +  l|I(x)  -  X(x*)|l]l|x  -  x*||+ 

+  ||i(x*)(x  -  X*)  +  Vf  (X-)  -  Vnx)Il} 

Moreover  from  (2.25)  and  proposition  3.2.5  p.  70  of  [19]  we  have 

(2.34)  ||9(x)||  <  a,||x  -  x*|l  V  x  6  S 

Hence  from  (2.30),  (2.25),  (2.33)  and  (2.34)  for  some  constants  aj,Q3  >  0 
we  have 

(2.35)  w(x,h)  <  otji  +  Oallx- x'll  VxgS,  Vh>h 

From  (2.27),  (2.31),  (2.32)  for  Xn,x„>i  6  5  and  >  h  we  have 


||x„+i-x*||  < 


<llG(x,,,h„)-x*||  +  ^j^||A;‘[(x„-x*)+(**-x,,_i)]y 

(2.36) 

< 

[(..(x„A)+^^^_J||x,  X  II  +  Jlx,-.  x|l< 

< 

«,e  +  «,  j  +  ^1  ll>t.  -  *•11  +  -  *'ll 

Moreover  from  (2.36)  eventually  changing  the  values  of  S  and  h  we 
have 

(2.37) 

.  O]  uo  1 

73  =  a3«  +  y+-^<  2 

ua  1 

^^  =  F<2 

so  that 


(2.38)  ||x„+i  -  x*i|  <  Tallxn  -  x*]]  +  74||Xn-i  -  x*|| 
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with  Os  =  73  +  74  <  1  that  is  Xn+i  6  S.  In  particular  we  have  shown  that 
(2.39)  lim  x„  =  x* 

fl^OO 


that  is  X*  is  a  point  of  attraction  of  (2.27). 

In  particular  for  n  >  Uq  >  0,  Xn  €  S,  using  (2.36)  the  required 
order-of- convergence  estimates  follows  from: 


(2.40) 


llXn+i  -  x*ll  <  +  OallXn  -  x*ll  llx„  -  x*l|  + 

+  T-T— llXo  -  x„_i||  for  n  >  no  >  0 

^n^n-1 


and  the  fact  that 


(2.41)  l|V/(jt„)ll  <  (lli(x-)ll  +  £)l|x,  -  x-|| 
where  lim  £„  =  0. 

n-*oo 

Using  the  method  given  by  (2.21),  (2.17),  (2.18)  requires  the  solution 
of  the  linear  system  (2.21)  at  each  step.  Computing  the  exact  solution 
with  a  direct  method  such  as  Gaussian  elimination  is  very  expensive  when 
a  large  number  of  unknowns  is  involved  and  may  not  be  worthwhile  when 
Xk  is  far  from  x*.  In  this  case  it  seems  natural  to  solve  the  linear  system 
(2.21)  by  an  iterative  procedure  amd  to  accept  an  approximate  solution. 
In  particular  since  the  matrix  A„  is  symmetric  and  positive  definite  we 
may  use  conjugate  gradients.  When  the  method  given  by  (2.21),  (2.17), 
(2.18)  is  used  to  solve  (2.21)  with  an  iterative  procedure,  accepting  an  ap¬ 
proximate  solution,  we  will  describe  this  procedure  as  an  inexact  method. 

Let  s„  be  the  approximate  step  computed  by  the  iterative  procedure 
when  solving  (2.21)  and 

(2.42)  r„  =  AnS„  -  b„ 

be  the  residual.  When  r„  =  0  the  linear  system  is  solved  exactly.  Let 
us  assume  that  the  approximate  step  computed  s„  satisfies  the  following 
condition: 


(2.43) 


l|r„||  <  0n\\K\\ 


n  =  0,1,... 
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for  some  forcing  sequence  {Pn},  n  =  0, 1,. . .  We  have  the  following  theo¬ 
rem: 

Theorem  2.2.  Let  f  :  IR^ — >01^  6c  twice  continuously  differ¬ 
entiable,  F{x)  =  f(x)^f(x)  and  L{x)  be  given  by  (2.15).  Let  x*  6  IB.^ 
6c  such  that  f(x*)  =  0,  ./(x*)  is  nonsingular  and  the  following  Lipschitz 
condition  holds: 

(2.44)  lli(x)  -  I(x*)||  <  7l|x  -  x*l|  V  X  e  S  =  {x|  llx  -  'c'll  <  6} 

for  some  constants  7,6  greater  than  zero.  In  the  iteration  (2.21),  (2.17), 
(2.18)  let  {h„},  n  =  0,1,2, ...,  6c  a  sequence  of  positive  numbers  and  let 
the  linear  system  (2.21)  6e  solved  approximately  in  such  a  way  that  the 
residuals  r„  given  by  (2.42)  satisfy  the  condition  (2.43)  for  some  forcing 
sequence  {0„},  n  =  0,1,....  If  0  <j3„  <  ,5max  <  1,  n  =  0, 1,...,  then 
there  exists  h  >  0  such  that  t/  /i„  >  6,  n  =  0, 1, . . .,  then  x’  is  a  point  of 
attraction  of  the  inexact  method  (2.21),  (2.17),  (2.18). 

Proof.  Since  7(x*)  is  nonsingular  and  X(x*)  =  2J{x‘)^ J{x‘)  we 
define  the  following  norm: 

(2.45)  11x11.  =  lll(x*)x|l  VxelR^ 
we  have 

(2.46)  ^IWI  <  W|.  <  miwi  VxeiR" 
where 

(2.47)  ;i,=max{lll(x*)ll,lll(xr‘ll} 

Moreover  it  is  easy  to  see  that  under  the  stated  hypotheses  for  any 
e  >  0  there  exists  6  >  0  and  h  >  0  such  that: 

(2.48)  ll>l(x,6)  -  X(x‘)ll  <  £  V  X  €  S  =  {x|  llx  -  x’H  <6},h>h 

(2.49)  ll44(x,  6)-»  -  I(x-)-‘  11  <  £  V  X  €  S  =  {x|  llx  -  x*ll  <  6}  ,  6  >  6 
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(2.50) 


||VF(x)  -  Vf (X-)  -  I(x*)(x  -  x*)|l  <  £||x  -  x*|| 


V  X  e  S  =  |xj  Ijx  -  x||*  < 


We  have 
(2.51) 

I(x*)(x„.,,  -  X*)  =  [/  +  I(x*)(yi;‘  -  I(x*)-i)]  • 

■{r„  +  (A„  -  I(x*))(x„  -  X*)  -  [  -  b„  -  VF(x*)  -  i(x*)(x„  -  x*)]} 


and  taking  norms: 


(2.52) 


l|in+,  -  x-||.  <[l  +  ||i(x-)|l||A;'  -  I(x-)-'l|]. 

■[ll'■n|l  +  P.  -  i(*')ll  ll*n  -  X'l|  + 

+||-b,-V/(x-)-i,(x-)(x,-x-)||] 


from  (2.24)  if  x^  €  5  and  /i„  >  h  using  (2.48),  (2.49),  2.50)  we  have: 


llx„+i  -  x-11.  <  [l  +  Mi^]  [4n||VF(x„)l|  +  2€llx„  -  x*|l+ 

(2.53) 

+  (1  +  4n)(||x„  -  X*||  +  llx*  -  X„-,||)] 

moreover  from 

(2.54)  VF(x,.)  =  i;(x*)(x„-x-)  +  [VF(x,.)- V/’(x*)-£(x*)(x,.-x*)] 
we  have 

(2.55)  ||VF(x„)||  <  |lx„  -  x*|l.  +  £||x„  -  x*l| . 

Finally  from  (2.47),  (2.53),  (2.55)  we  have: 

Il*"+1  -  ^m«(l  +  £/ii)  +  (2  +  ^  • 

•||X„  -  X*||.+[l  +/tlfl^(l  +  /3ma.)||x„„,  -  X*||  = 
=  QsllXn  -  X*||.  +  OfillXn.!  -  x‘|| 


(2.56) 
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where 


(2.57) 


Os  =  [1  +  /iif] 


0!6  =  (1  +Plf)(l  +  0tn»x 


choosing  the  values  of  £  and  h  so  that  aj  +  <  1  frorr  (2.56)  we  have 

that  if  x„,x„_i  6  S  then  x„+i  e  S  and 


•’m  x„  =  X* 
—  00 


Theo-REM  2  3.  let  f  ;  IR.^ — ►IR^  6e  twice  continuously  differ¬ 
entiable,  F{x)  =  f(x)’’f(x)  and  L(x)  be  given  by  (2.15).  Let  x*  € 
bi  such  that  {(^‘)  =  0,  J'(x’)  is  nonsingular  and  the  following  Lipschitz 
condition  holds. 

(2.58)  l|ii(x)  -  I(x*)|l  <  7l|x  -  x‘11  V  X  5  =  {xl  Hx  -  x*ll  <  S) 

In  the  iteration  (2.21),  (2.17),  (2.18)  let  {/»«}.  ^  =  0,1,...,  be  a 
sequence  of  positive  numbers  and  let  the  linear  system  (2.21)  he  solved 
approximately  in  such  a  way  that  residuals  r„  given  by  (2.42)  satisfy  the 
condition  (2.43)  for  some  forcing  sequence  {P„,  n  =  0,1,...},  such  that 
0  <  4«»  <  4m*x  <  1(  n  =  0, 1, —  Then  there  exists  h  such  that  if  hn  >  h, 
n  =  0,1,..., X*  is  a  point  of  attraction  of  the  inexact  method  (2.21), 
(2.17),  (2.18)  and  the  rate  of  convergence  is: 

(i)  Q-superlinear  if  h~^  <  7i||VF(x„)||,  7i  >  0,  n  >  no  for  some 
7i,no  >  0  and  lim  /?„  =  0 

n-*co 

(ii)  Q -quadratic  if  h'^  <  7jl|VF(Xn)ll*,  7,  >  0,  n  >  no  and  4n  < 
7allVF(x„)ll,  7j  >  0,  n  >  no,  for  some  7,,  no  >  0. 
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Proof.  From  Theorem  2.2  we  have  that  x*  is  a  point  of  attraction 
of  the  inexact  method  (2.21),  (2.17),  (2.18)  so  that  we  can  assume  that 
lim  x„  =  X*  and  it  remains  to  prove  the  rate-of-convergence  results. 

fl  >00 

We  have: 


^+1  -  X*  =  -  Z(x*)|(x„  -  x*)+ 

(2.59) 

-[-b„-VF(x*)-I(x*)(x„-x*)]} 

and  taking  norms 

lli.+i  -=c-|l  <  M;'ll[llr.||  +  IM.  -  i(x-)||  ||x.  -X-II+ 
(2-60)  +  llvnx.)  -  Vf(x-)  -  I(x-)(x„  -  x-)||+ 


Let  £,6,h  be  chosen  in  such  a  way  that  (2.29),  (2.34),  (2.35)  hold 
then  there  exists  Uq  such  that  for  n  >  nj,  +  1,  x„  6  S  =  {x]  jjx  -x*|l  < 
we  have: 


(2.61) 


||x„+i-x*|l<a[4„||VF(x,.)ll+ 

+  (<»2)^+“3llx„  -x’ll^llXn  -X’||  + 


+  “lliXn  -  X*|i^  +  ^(1  +  /j„)(|Xn  -  X„_i|I 
and  the  desired  rate-of-convergence  results  follow  from  (2.41). 


3  -  Complementarity  problems  and  nonlinear  systems 

Let  f  :  IR^  — ►  be  given,  the  complementarity  problem  associ¬ 

ated  with  f  is 

(3.1)  X  >  0 

(3.2)  f(x)  >  0 

(3.3)  {x,f(x))  =  0 
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and  let  0  :  IR  — *■  IR  be  a  strictly  increasing  function  such  that  0(0)  =  0. 
In  [17]  Mangasarian  has  shown  that  x*  G  IR^  is  a  solution  of  the 
complementarity  problem  (3.1),  (3.2),  (3.3)  if  and  only  if  x*  is  a  sdution 
of  the  system  of  nonlinear  equations 

(3.4)  g(x)  =  0 
where  g(x)  =  (5i(x), ja(x),...,5A,(x))^  and 

(3.5)  gi{x)  =  0{\fi(x)  -  Xi\)  -  0{fi(x))  -  9{xi)  i=  1,2,... ,1V 
for  later  purposes  let  us  introduce 

(3.6)  G(x)  =  g(x)’-g(x) 

Definition  3.1:  Let  x*  €  IR^  be  a  solution  of  the  complementarity 
problem  (3.1),  (3.2),  (3.3)  we  will  say  that  x*  is  nondegenerate  if  x*  + 
f(x*)  >  0. 

Definition  3.2:  Let  f  be  continuously  differentiable  and  J{x)  = 
df  /dx  be  the  Jacobian  of  f  with  respect  to  x,  if  for  n  =  1, 2, . . . ,  IV  each 
principal  minor  {{dfildxj)),i,j  =  1,2, ...,n,  is  nonsingular  we  say  that 
7(x)  has  nonsingtilar  principal  minors. 

In  [17]  Mangasarian  has  shown  that  if  x*  is  a  nondegenerate  solu¬ 
tion  of  the  complementarity  problem  (3.1),  (3.2),  (3.3)  such  that  J{x') 
has  nonsingular  principal  minors  and  0  :  IR  — ►  IR  is  a  strictly  increasing 
differentiable  function  such  that  d0/dt{O)  -|-  d0/dt{t)  >  0,  Vt  >  0,  then  x* 
is  a  solution  of  the  nonlinear  system  (3.4)  and  5g/5x(x*)  the  Jacobian  of 
g  with  respect  to  x  is  nonsingular. 

For  simplicity  we  choose  0{t)  =  t/2  so  that  in  a  neighborhood  of 
a  nondegenerate  solution  of  the  complementarity  problem  (3.1),  (3.2), 
(3.3)  the  function  g(x)  given  by  (3.5)  has  the  same  regularity  properties 
of  f(x).  Given  the  local  character  of  the  convergence  theorems  of  section 
2  this  is  satisfactory.  In  section  4  the  method  for  solving  nonlinear  system 
described  in  section  2  will  be  applied  to  (3.4)  with  0{t)  =  t/2  for  some 
test  complementarity  problems. 
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4-  Numerical  experience 

The  inexact  method  (2.21),  (2.17),  (2.18)  has  been  implemented  as 
follows: 

i)  since  An  is  symmetric  and  ptwitive  definite  the  linear  system  (2.21) 
has  been  solved  by  the  conjugate  gradient  method  (C.G.)  introduced 
by  Fletcher  and  Reeves  [20].  This  procedure  solves  aji  N  x  N 
linear  system  in  at  most  N  steps.  Hovewer  we  stop  the  conjugate 
gradient  procedure  after  a  number  of  steps  which  is  usually  consider¬ 
ably  lower  than  N.  In  fact  let  be  the  approximate  value  for  the 
solution  s„  of  the  linear  system  (2.21)  obtained  as  the  result  of  step 
k  of  the  conjugate  gradient  iteration;  the  iteration  is  stopped  after 
step  m  if 

IKs(„"‘)-b„l|</3„||b„|| 

ii)  We  have  chosen: 

^0  =  »7o  =  0 
A*  =  5  =  1 

D  =  I  (the  identity  matrix) 
s^°)  =  0  n  =  0,l,... 

and  the  following  very  simple  variation  laws  for  the  time  integration 
step-length  h„  and  the  forcing  sequence  Pn- 

hn  +  l  =  niin(10/ln>  ^m»x)  71  =  0,  1, 2, .  .  . 

with  ho  =  1,  Amax  =  10” 

PU,  =  an0l  n  =  0,1,2,... 

where  /3o  is  given  and  a„  is  automatically  chosen  by  the  program 
among  the  two  values  0.1  and  0.5. 

iii)  The  program  stops  in  any  case  the  conjugate-gradients  iteration  after 
A'  steps  in  order  to  avoid  possible  non  termination  due  to  the  finite 
arithmetic  of  the  computer. 
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Finally  the  method  given  by  (2.21),  (2.17),  (2.18)  (i.e.  exact  solution 
of  the  linear  system  (2.21))  is  obtained  simply  setting  jdo  =  0. 

The  stopping  rule  adopted  is  G'(x„)  <  10“*°  for  the  inexact  method 
and  G(x„)  <  10“*“  for  the  “exact”  method  (i.e.  /So  =  0).  These  methods 
have  been  coded  in  the  Pascal  programming  language  and  the  program 
has  been  run  on  a  Hewlett-Packard  9816  computer. 

We  have  tested  the  proposed  algorithm  on  three  complementarity 
problems  of  which  two  are  linear  and  one  is  nonlinear. 

The  first  problem  considered  arises  as  a  one-dimensional  free-boun- 
dary  problem  in  the  lubrication  theory  of  an  infinite  journal  bearing, 
i.e.  a  rotating  cylinder  separated  from  a  bearing  surface  by  a  thin  film  of 
lubricating  fluid  [21].  The  finite-difference  approximation  used  by  Cryer 
in  [21]  leads  to 

Problem  A  (called  Problem  3D  by  Cryer):  Find  x,  w  e  IR^  such 

that 


(4.1)  w  =  q+Mx,  w>0,  x>0, 

(4.2)  (w,x)  =  0 

where  M  =  1, 2, . . . ,  JV,  is  an  x  JV  matrix  with  elements 

Mij  given  by 


(4.3) 


Af.-;  =  -(H.>:/2)^ 

if  j  =  i+l. 

M,,  =((H.>x/2)^-H(H._i/2)1, 

11 

Mx^  =  -(H._x/,)“, 

if  ;■  =  i  -  1 , 

o 

II 

otherwise 

and  q  =  (?!, ?2» •  •  • » 9w)^  is  a  vector  with  elements  given  by 


(4.4) 
where 

(4.5) 


T 

9i  =  j[^i+l/2  -H,-i/2],  »  =  1,2,. ..iV 


•^i±l/2  +  B 


T 

iV  +  1. 
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and  the  function  H{y)  is  given  by 

(4.6)  E{y)  =  -^{1  +  e  cos  ry)  >  0 
with 

(4.7)  T  =  2,  £  =  0.8 

We  note  that  the  matrix  M  given  by  (4.3)  is  symmetric  and  positive- 
definite. 

The  second  problem  arises  as  a  two-dimensional  free-boundary  prob¬ 
lem  in  the  theory  of  the  steady-state  fluid  flow  through  porous  media. 
Some  of  these  problems  can  be  formulated  as  a  variational  inequality  af¬ 
ter  an  ingenious  transformation  proposed  by  Baiocchi  and  others  (ref. 
[13]).  The  discretization  used  on  the  “model  problem”  ([3],  p.  4)  leads  to 

Problem  B:  Find  x,  w  6  such  that 

(4.8)  w  =  q-f-A/x,  w>0,  x>0, 

(4.9)  (w,  x)  =  0 

where  M,  an  iV  x  TV  real  matrix,  and  q  =  (91,92, ... ,g;v)’^  G  IR^  are 
defined  below. 

Given  n,,n,  (positive  integers)  and  X,  Y  (positive  real  numbers),  let 

N  =  rifTiy , 

Dx  =  Jf/n,  +  1 , 

Dy  =  Y/riy  +  1 , 
a  =  DyJDx , 

let  A  be  the  n,  x  n,  tridiagonal  matrix  having  all  the  main  diagonal 
elements  equal  to  2(a  +  1/a),  and  the  paradiagonal  elements  (i.e.  im¬ 
mediately  above  or  below  the  main  diagonal)  equal  to  -a,  and  let  B  be 
the  Tig  X  rig  diagonal  matrix  with  diagonal  elements  equal  to  —  1/a.  The 
matrix  M  is  an  n  x  n  matrix  with  a  block-tridiagonal  structure  (n,  x  n. 
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blocks),  having  each  main-diagonal  block  equal  to  the  matrix  A,  and  each 
paradiagonal  block  equal  to  the  matrix  B.  We  note  that  A/  is  a  positive- 
definite  symmetric  matrix.  The  vector  q  is  defined  as  follows.  Given  W 
(0  <  W  <  y),  and  using  the  Kronecker  symbol  let 

u(i)  =  1(1'  -  if . 

9M=l(W-y)\  if  y<W, 

9R{y)  =  0 ,  if  y  >  w , 

PD(x)  =  yV2-(y^-w^)(x/2jr), 

gu{x)  =  0, 

Tij  =  -DxDy  +  6iiagi,{jDy)  +  Si„,agR{jDy)+ 

+  6ij{l/a)gD{iDx)  +  6n,jilfa)guiiDx) , 

i  —  1,2,..., Tij, ,  j  ^  1, 2, . . . , Tij  . 

The  elements  gi ,  92»  •  •  • » 9n  of  q  are  given  by 

(4.10)  qk  =  Tij ,  with  k  =  {j  -  l)n,  -|-  i 

Our  last  problem,  which  is  defined  below,  can  be  interpreted  as  a 
finite-difference  approximation  of  a  nonlinear  variational  inequality. 

Problem  C:  Find  x,w  e  such  that 

(4.15)  w  =  Mx  4- p(x) -i- q ,  w>0,  x>0 

(4.16)  (w,x)  =  0 

The  problem  dimension  N,  the  quantities  Dx,  Dy  and  the  matrix 
M  are  defined  as  in  problem  B,  given  n,,n,,A’,y.  The  nonlinear  term 
p(x)  is  a  vector  in  IR^  with  components  p,  =  xf,  i  =  The 

vector  q  =  (qi,9j, . . .  ,?n)^  is  defined  by  equation  (4.10)  where  = 
DxDy  s\a(2tiDx / X),  i  =  1,2,. . . ,n*,  j  =  1,2, ...  ,n,. 

The  numerical  results  obtained  with  the  previously  described  meth¬ 
ods  on  Problem  A,  B,  C  are  shown  in  Table  1,  2,  3  respectively. 


[21] 


An  Inexact  Continuous  Method  for  the  Solution  etc. 


541 


542 


P.  ALUFFI-PENTINI  -  V.  PARISl  -  F.  ZIRILU 


{22) 


In  tables  1,  2,  3  the  adavantageof  using  “inexact  Unear  algebra”  with 
respect  to  complete  solution  of  the  Unear  system  for  problems  A,  B,  C  is 
shpwn,  and  the  advantage  is  increasing  with  the  number  of  unknowns. 
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abstract 

.A  new  method  to  solve  linear  programming  problems  is  introduced.  This  method 
follows  a  path  defined  by  a  system  of  o.d.e..  and  for  nunJegenerate  problems  is 
quadratically  convergent. 


1.  INTRODLCTION 

Let  R"  be  the  n-dimensional  real  Euclidean  space,  and 

x=(jCi . x„)^^R".  where  the  superscript  T  means  transpose.  For  x, 

y  e  fl"  let  x’^y  be  the  usual  Euclidean  inner  product,  and  let  e  =(1 . 1)^ 
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A  linear  programming  problem  consists  in  minimizing  a  linear  ionction  oiver 


a  region  defined  by  linear  equality  and  inequality  constraints.  We  will  say 
that  a  linear  programming  problem  is  in  canonical  form  when  it  is  written  as 
follows: 

minimize  c^x 
»e  H- 

(1.1) 

subject  to 

Ax  =  0, 

(1.2) 

II 

(13) 

X  >  0, 

(’.4) 

with  side  conditions  Ae  =  0,  where 

A  =  ((fl,,))  1  .  ,  .  n  >  2, 

m  <  n. 

and  c  €  R"  are  given,  and  the  inequality  (l.-l)  is  understood  componentwise, 
that  is, 

ij^O,  _/  =  1 . n . 

Moreover,  we  assume  that  the  matri.x  A  is  of  rank  m. 

We  note  that  on  these  hypotheses  (l/n)e  is  a  feasible  point,  so  that  the 
feasible  region  is  not  empty  . 

The  simplex  method  applies  to  linear  programming  problems  in  standard 
form,  that  is. 


minimize  d^y 

(15) 

y  €  H' 

subject  to 


Cy  >  b. 

(1.6) 

y  >  0. 

(1.7) 

where  d  e  fi’,  b  e  fl',  C  e  fi'’*’.  r  <  s.  are  given,  and  the  inequalities  (1.6), 
(1.7)  are  understood  componen^vise.  In  (l)  it  has  been  shown  that  a  linear 
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programming  problem  in  standard  form  with  a  finite  solution  ean  alwaxs  be 
reduced  to  canonical  form.  Moreoxer.  fCannarkar  assumes  that  the  objectixe 
function  of  the  problem  (1  1)-(1.4)  is  such  that 

z»  =  cV=0 

for  anx-  feasible  point  x*  that  is  a  solution  of  the  linear  procramming  problem 
The  problem  (1.1)-(1.4)  xxith  this  extra  assumption  is  called  a 
problem  in  canonical  form  with  a  normalized  objective  function.  In  ou.'  xvork 
The  assumption  of  haxing  a  normalized  objectixe  furction  is  not  neeessarx, 
ho  *ever,  since  this  assumption  simplifies  some  of  the  following  algebraic 
manipulations,  xve  xxill  keep  it. 

Let  n  denote  the  subspace  0  =  (x  e  R“  |  .Ax  =  0).  let  A  lx*  the  simplex 
A  =  {x  G  fi"  lx  >  0.  e^x  =  U.  and  finally  let 

.\  =  anA  (l.S) 

be  the  polytope  of  the  feasible  points.  Then  the  problem  (1.1)-(1.4)  can  be 
rexx  ritten  us  follov  s: 


minimize  c^x. 
« e  \ 


(1.9) 


In  this  paper  xve  xxill  introduce  a  nexv  method  to  sobe  linear  program¬ 
ming  problems  in  canonical  form  xxith  a  normalized  objective  function.  This 
class  of  problems  is  the  one  considered  by  N.  Kamiarkar  in  his  celebrated 
iwper  [2]. 

In  the  late  1940s  C.  B.  Dantzig  [3j  developed  the  simplex  method  to 
solve  linear  programming  problems.  In  1972  V.  Klee  and  C.  L.  .Mints  [4] 
shoxved  that  the  worst  case  complexity  of  the  simplex  method  is  combinato¬ 
rial.  Here  the  term  ■‘complexicx  ”  means  the  number  of  elementary  operations 
necessary  to  soixe  a  linear  programming  pn>blem  in  the  standard  Ibrm 
(1.5)-(1.7).  Sinc-e  the  simplex  method  finds  the  solution  after  a  finite  numlier 
of  iterations.  Klee  and  Minty  (4]  xvere  able  to  give  an  example  where  ..ie 
simplex  method  has  complexity 


p  =  0{rs2'). 


Note  that  in  (1.5)  y  6  R'.  .Moreover,  in  1981  S.  Smale  in  [5]  showed  that  the 
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"average  ”  complexity  of  the  simplex  method  is 

p  =  0(n*). 

where  r,s  are  the  dimensions  of  the  matrix  C  in  (1.6). 

In  spite  of  its  worst  case  combinatorial  complexity,  the  simplex  method 
has  been  very  successftjl  in  solving  linear  programming  problems.  The 
feature  of  the  simplex  method  responsible  for  its  worst  case  combinatorial 
complexiri  is  that  it  moves  on  the  boundary  of  the  feasible  region 


Q  =  {y  €  fi’ICy  >  b.  y  >  0} . 


In  recent  years  a  great  deal  of  effort  has  been  spent  in  the  attempt  to  find  a 
new  algorithm  for  linear  programming  whose  coniple.xih  in  the  worst  case  is 
polynomial.  It  is  believed  that  these  new  methods  will  go  through  the 
interior  of  the  feasible  region  Q. 

In  1979  L  C.  Khachijan  f6j  introduced  the  first  method  of  this  class, 
called  the  ellipsoid  medxod.  The  worst  case  complexity  of  this  method  is 

P  =  0(s") 

Here,  however,  the  meaning  of  the  term  "eomplexitv"  has  been  slightly 
changed.  In  fact  the  ellipsoid  method  does  not  step  after  a  finite  number  of 
iterations,  so  that  "complexity  means  the  number  of  elementarv  operations 
neeessarv-  to  arrive  in  a  predetennined  neighborhood  of  the  solution.  More¬ 
over,  the  method  introduced  by  Khachijan  is  only  of  theoretical  interest, 
since  its  practical  perfonnance  is  rather  poor. 

In  I98-I  .\.  Kaniiarkar  [2j  presented  a  new  linear  programming  method  of 
polynomial  worst  case  complexity 


P  =  0(s”). 

This  algorithm  is  called  the  projective  method  when  applied  to  a  linear 
programming  problem  in  canonical  form  with  a  normalized  objective  func¬ 
tion.  This  algorithm  is  of  theoretical  and  practical  importance. 

Since  1984  a  great  deal  of  work  has  been  done  in  developing  new 
methods  for  linear  programming.  Several  “interior  point  algorithms'  have 
been  proposed.  P.  E.  Gill.  W.  .Murray.  M.  .V  Saunders.  J.  A.  Tomlin,  and 
M.  H.  Wright  in  (7]  have  interpreted  Karmarkar  s  algorithm  as  a  "logarithmic 
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barrier  method”  an  I  have  suggested  a  new  algorithm  with  good  practical 
performance. 

In  the  framework  of  logarithmic  Iwrricr  function  methods  we  can  recall 
the  work  of  several  authors.  In  (8)  J.  Renegar  lowered  Karmarkar's  ct)m{^- 
it\’  bound.  In  (9]  C.  Conzaga  lowered  Renegar's  comple.vitv'  bound.  In  (lOj  M. 
Iri  and  H.  Imai,  with  the  hspothesis  of  being  able  to  perform  exact  line- 
searches,  introduced  a  (juadratically  convergent  algorithm  fo.  the  linear 
programming  problem.  In  [ll]  N.  Megiddo  studied  the  geometrica.  pn.per- 
ties  of  the  paths  derived  from  “weighted  logarithmic  barrier  functions.” 
Finally,  j.  A.  Tomlin  in  (12)  reports  on  considerable  numerical  experimenta¬ 
tion  with  this  kind  of  algorithms. 

In  this  paper,  as  suggested  by  D.  .A.  Bayer  and  J.  C.  Lagarias  in  [13],  we 
will  show  that  Karmarkar’s  projective  method  can  be  obtained  b>-  appKing 
Euler  s  method  w  ith  variable  stepsize  to  a  suitable  initial  value  problenr  for  a 
system  of  ordinary  differential  erjuations.  In  fac-t-  Karmarkar  s  method  obtains 
the  solution  x*  of  the  linear  programming  proldem  by  computing 


where  x(f.(l/ rj)e)  is  the  solution  of  a  system  of  ordinary  differential  equa¬ 
tions  with  initial  condition  (l/n)e,  using  Euler’s  method  with  variable 
stepsize.  The  idea  of  obtaining  the  solution  of  nonlinear  programming 
problems  as  limit  points  of  the  trajectories  of  systems  of  ordinary  differential 
equations  has  been  widely  used;  for  a  review  see  [14j.  In  particular,  in 
[15-17]  quadratically  convergent  algorithms  for  nonlinear  systems  of  equa¬ 
tions  have  been  obtained  from  methods  based  on  the  numerical  integration  of 
trajectories  of  systems  of  ordinary  differential  equations. 

The  interpretation  of  Karmarkar’s  projective  method  as  the  numerical 
solution  of  an  initial  value  problem  raises  two  natural  questions: 

(i)  Can  the  system  of  ordinary  differential  equations  used  in  Karmarkar  s 
projective  method  be  changed  to  a  new  one  that  will  generate  an  interesting 
algorithm’:* 

(ii)  Can  the  Euler  method  with  variable  stepsize  that  is  used  in  Kar- 
markar's  projective  method  be  replaced  with  some  other  numerical  scheme 
that  will  generate  interesting  algorithms’ 

An  answer  to  question  (i)  has  been  given  by  D.  .A.  Bayer  and  J.  C. 
Lagarias  in  [13]  and  J.  L  Nazareth  in  [18],  who  replaced  Karmarkar’s  vector 
field  with  the  ajftne  vector  field.  Question  (ii)  has  been  considered  by  .N. 
Karmarkar.  J.  C.  Lagarias.  L.  Slutsman.  and  P.  W  ang  in  [19],  where  they  tried 
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to  approximate  the  path  x(t.(l/n)e)  with  a  power  series  expansion,  obtaining 
encouraging  practical  results.  In  this  paper  we  give  two  new  answers  to 
questions  (i)  and  (ii);  in  fact,  we  propose  a  vector  field  which  is  difierent 
from  the  ones  previously  considered,  and  we  use  a  linearly  implicit  A-stable 
integration  scheme  [14]  to  solve  the  initial  value  problem  considered.  In  this 
way  we  obtain  a  quadratically  convergent  algorithm  for  linear  programming 
problem.  .Moreover  our  algorithm  shows  good  practical  behavior. 

In  Section  2  Karmarkar's  projective  method  is  interpreted  as  the  numeri¬ 
cal  integration  of  an  initial  value  problem  with  Euler's  method  and  variable 
stepsize.  Moreoser.  to  a  linear  programming  problem  in  canonical  form  with 
normalized  objective  function  is  associated  a  new  system  of  ordinary  differ¬ 
ential  equations.  If  we  assume  that  the  solution  of  the  linear  programming 
problem  is  unicjue.  this  solution  can  be  obtained  as  the  limit  point  of  a 
suitable  trajectory  of  the  system  of  ordinary  differential  equations. 

In  Section  3  an  initial  value  problem  for  this  system  of  ordinary  differen¬ 
tial  ecjuations  is  integrated  numerically,  using  a  linearly  implicit  .A-stable 
method  with  variable  stepsize.  It  is  shown  that  this  is  a  (juadratieally 
consergent  algorithm  for  linear  programming. 

Finally,  in  Section  4  we  compare  the  computational  cost  of  our  step  with 
that  of  Karmarkar  s  projective  algorithm  and  that  of  the  simplex  algorithm, 
and  we  present  some  numerical  experiments. 

2.  THE  LSE  OF  ORDIN.ARY  DIFFERENTIAL  EQL'.ATIONS 
IN  LI.\E.AR  PROCaAMMINC 

Let  X*  =(i*.xr . x,:)’'e  fl".  A*e  fl""’  Ise  given  by  .V*  =  ((.V* ))  = 

Diaglx*),  that  is.  .V *  =  x*5,^,  i.j  =  1,2 . n.  where  5,^  is  the  Kronetker 

symlx)!. 

DtKiMTiov  2.1.  .A  minimizer  x*  of  the  problem  (1.1)-(14)  is  called 
nondegeiierate  if  it  has  exactly  «  —  wi  —  I  null  ecmiponents. 

bet  7„={1.2 . «}  and  S={s,.5. . s„,, ,).  m 1  ^  n.  be  an  ordered 

set  of  indices  such  that  S  C  y„.  Let  z  =(z,.;, . -,)€  R"  be  a  vector.  We 

denote  by  z.,  the  sector  z.,  =(r,  . ^)e  R'“~'.  .Moreover,  gisen  a 

vector  V  e  R"'~'  and  a  matri.x  ^  denote  by 

the  submatrix  =[</’'.</'- . «/’■•*')€  flo-i- i>»(m- u  ij  jjp, 

column  of  Q.  If  B  is  an  ordered  set  of  indices  such  that  Bc]^  and 
.V  =  }„-  B.  then  the  system 


^z  »  V 


(2.1) 
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can  be  rewritten  in  the  following  form: 

+QVZV  =  '  (--) 

Definition  2.2.  Let  be  an  ordered  set  of  m  +  1  indices.  Then  B  is  a 
set  of  basic  indices  for  the  system  (2.2)  if  there  e.\ists  a  matri\  e 
^(in.>.i)x(n-fn-ii  ^  vcctor  V  G  R“~*  such  that  the  system  (2.2)  is  equiva¬ 
lent  to  the  system 


^8  ^ 


(2.3) 


Lemm.a  2.3.  Let  B  he  an  ordered  set  of  m  +  1  indices  such  that  B  is  a  set 
of  basic  indices  for  the  system  (2.2).  Then  Qg  is  an  invertible  matrix. 


Proof.  It  follows  immediatelv  from  the  eijuixalence  of  the  linear  s\  stems 

(2.2)  and  (2.3).  ■ 

Let  a>0.  x>0.  x“=(.t;.i? . i“)^efl\  and  be  the 

matrix  A'“  =  Diagix”). 

Lemm.x  2.4.  Let  o>l  and  x*  be  a  nondegenerate  minimizer  of  the 
linear  programming  problem  (l.I)-(l.4).  Then  is  an  invertible  ma¬ 

trix. 


Proof.  Let 


A/  = 


be  the  matrix  A  with  the  extra  row  added.  Since  x*  is  a  nondegenerafe 
minimizer  of  the  problem  (1.1)-(1.4)  and  .Vf  has  rank  m  +  1,  then  there 
exists  an  ordered  set  of  m  +  1  indices  B  such  that  x  J  e  fl  '  has  all  nonzero 
components  and  x  *  €  B""'""'  is  the  zero  vector,  where  .V  =  J^  —  B.  More¬ 
over.  B  is  a  set  of  basic  indices  for  the  system 


.Wy  -  u. 


(2.4) 


where  u  e  fl'"*' 


is  given  and  y  e  ft".  Let  ,v  e  i> 
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matrix  =  Diagixj'^®).  and  x*'/®  6  "  be  the  matrix 

Xy'^’®  =*  Diag(x*‘'^®),  that  is.  the  null  matrix.  So  from  Lemma  2.3,  AfjS 
^<m*i)x(«-n  jj  invertible,  which  implies  that  A/gX^*^®  is  invertible.  More¬ 
over,  since  B  is  a  set  of  basic  indices  for  the  system  (2.4),  we  have 

MX*.Vf®‘=  (.WaXg*'^®  +  .V/vX?'''®)(x;'^®.VfJ  + 

=  (.Vfg.Y;'/®)(Xg-/®.M;).  (2.5) 

Since  .VfgXg ''^®  is  invertible,  from  (2.5)  it  follows  that  is  invertible, 

so  that  an  easy  computation  shows  that  .AX*.A®^  is  invertible.  Therefore  it 
follows  that  AX**''®  is  of  rank  m.  so  that  .AX*"''®  is  of  rank  m  and  .A.Y*“A®^ 
is  invertible.  ® 


Lemma  2.5.  Let  a  =  1  or  a  =  2,  and  let  x*  be  a  nondegenerate 
minimizer  of  the  linear  programming  problem  (1.1)-(1.4).  Then  there  exists 
p*  >  0  such  that  AX‘‘A^  is  invertible  for  x€S(x*,p*),  where  S(x*,p*)  = 
{xefl-||lx-x*||<p*). 

Proof.  The  proof  follows  immediately  from  the  continuity  of  .A-V^A*^ 
with  resjject  to  x  €  R".  from  Lemma  2.4.  and  from  j.  M.  Ortega  and  U.  C. 
Rheinboldt  (20.  Lemma  2.3.2.  p.  45).  ■ 

Let  X  =(x,..t, . e  fi";  let  X  e  fi"*"  be  given  by  X  =  Diag(.x);  let 

D  €  fi"*",  m  <  n.  be  a  matrix,  and  D  *  the  subspace 


D- ={xefl'’|Dx  =  0};  (2.6) 

and  let  flo-t  )  be  the  orthogonal  projection  on  D~  ■  The  projector  flo-l  ) 
always  exists,  and  if  D  has  full  rank  is  given  by 

Ho-fy)  =  (/- D'’(DD'’)‘'D]y.  y  e  fl".  (2.7) 


Let  r  be 


r  =  {xefl'’|x>0). 


(2.8) 
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and  f  be 


its  interior: 


r  =  {xefl"U>0).  (2.9) 

Tbe  set  F  is  called  the  positive  orthant. 

For  a  >  0  and  x  e  F  let  X"'"  e  fl"''"  be  the  mafri.v  X"  ’  = 
Diag(if''^,iJ''^,....i“'-).  We  observe  that  n,^^.  can  ahvavs  be  e.\- 
pressed  in  the  form  (2.7).  In  fact,  let  r  be  the  rank  of  the  matri.v  AX“  if 
r  -  m,  Aen  11,^;^..,:,*  is  given  by  (2.7).  If  0  <  r  <  m.  we  can  consider  the 
matrix  A  €  obtained  from  A  by  eliminating  the  in  —  r  rows  of  .A  with 
indices  equal  to  those  of  the  _^m-r  rows  of  that  are  linearly 

dependent.  Since  (.A,V“'*)*  =(.AX“'")'.  we  have 

^(.AV*  ^1-  (2.10) 

where  fl, , is  given  by  (2.7).  Finally  if  r  =  0  we  have  that  11,  ^^..  =  /, 

where  /  is  the  n  X  n  identitv  matrix.  Let  h(x)€  fi"  be  the  follow  ine  vector 
field: 

h(x)  =  -.V(f-ee’‘X)n,.,.v,-(Xc).  x6fl".  (2.11) 

The  vector  field  h(x)  is  known  as  Karmarkar’s  vector  field  [2.  13).  Let 
S(x*,p*)  be  the  ofien  ball  of  Lemma  2''  and  let  us  consider  h{x)  for 
*  s  F  U  S(x*,p*).  ^ 

We  observe  that  for  x  s  FuS(x*,p*).  AX  is  of  rank  m.  Then  h(x)  is  a 
continuously  djfferentiable  function  of  x  for  x  6  Fu  S(x*,p*).  Let  \  be  given 
by  (1.8),  and  .\  be  given  by 


A  =  Anf.  (2.12) 

W'e  will  consider  the  initial  value  problem 


dx 

(2.13) 

x(0)  =  -e. 

(2.14) 

n 


It  is  easy  to  verify  that  (l/n)e  e  A. 
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Lemma  2.6.  Let 


B  = 


eB‘'" 


l)XJt 


be  the  matrix  AX  with  the  extra  row  added,  and  let  y  =  */vn(  n  - 1) . 
where  x  €  (0. 1)  is  a  parameter.  For  €  A  let 


Xi  =  Diag(x*). 


and  let  be  given  by 

A<*»|;;^i|n8^.(.V)||-e%n3,.(XiC)j  .  (2,15) 

Then  Euler's  method  applied  to  the  initial  value  problem  (2.13)  12.14)  with 
variable  stepsize  Alj  given  by  (2.1.5)  produces  the  sequence  (x‘).  k  = 

0.1,2 .  generated  by  Karmarkar's  algorithm  (2.  pp.  37S-379j  applied  to 

the  linear  programming  problem  with  normalized  objective  function 
(1.1)-(1.4). 


Proof.  We  observe  that  Atj  >0  for  e  .\  (see  [2.  Theorem  5.  pp 
381—382],  so  that,  integrating  (2. 13).  (2.14)  with  Euler’s  method  and  variable 
stepsize  A/j.  we  have 


II 


(2  16) 


x‘*'  =  x‘  +  Atth(x‘)  it=0.1.2. 

The  thesis  follows  from  a  straitshtfoiward  computation 
For  a  >  0  let 


(2  17) 


g(x.a)  =  n,,v*  :,.( .V'  -c). 


xe  r 


(2  18) 
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We  note  that  for  a  =  2,  g(x.2)  is  the  affine  scaling  factor  of  [13],  and 

g(x.2)  =  n,,.v,-(-Vc)  (2.19) 

can  be  defined  for  x  e  fi"  so  that 

h(x)  = -A•(/-ee^V)g{x.2).  x  e  fiV  (2  20) 

We  have: 


TiiEOBtM  2.7.  For  0  <  o  <  2  /cf  x  6  A  be  a  feasible  point  for  the  linear 
programming  problem  with  normalized  objective  function  (1.1)-(1.4).  Then 

n,,v.  m-(A“  'c)  =  0  (2  21) 

if  and  only  if  x  is  a  tninimizer  of  the  linear  protiramniinf’  problem  with 
normalized  objective  function  (1.1)~(1.4). 


Proof.  Let  x  =  Xe  be  a  minimizer  of  the  problem  (l.l)-(l  4)  with 
normalized  objective  function.  Then 

(2.22) 

In  fact,  if  we  assume  that 

n,.*.v.:,-(A"^'c)#0.  (2.23) 


then  there  exists  z  =(Z|,Zj . -,)€  R"  such  that 

53  -0.  i=l . m.  (2.24) 

that  is.  z  ^  and 


c 


p*o. 


(2.25) 
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that  is.  is  not  orthogonal  to  z.  We  can  assume  without  loss  of 


generality  ^  >  0.  Let  us  define 

Wj~x“^^Zj.  j-1.2 . n.  (226) 

Since  0,  there  exists  j  such  that  itv  #  0.  If  Wj  <  0  for  j  =  1.2 . n.  we 

choose  e  >  0;  otherwise  we  choose  e  as  follows: 

0<c<  min  — .  {227) 

j  ttj  >  0 


We  recall  that  >  0  for  =  1.2 . n.  Let 

t-;  =  ;  =  1,2 . n.  (2  28) 


From  (2.27)  we  have 


Vj>Q.  j  =  1.2 . n. 


and 

ft-  >0. 

Let  us  define 


(2.29) 


(230) 


r 


1.2. 


I 


(2.31) 


The  point  u  “(uj.u^ . *  feasible  point;  in  fact, 


•Au  -  0. 

(232) 

e’^u  =  1. 

(2  33) 

u  >  0 

(2  34) 
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Moreover. 


«  1  * 

^  ^j'‘j  “  ~  E 

j-i  ®  ''  ;-i 

=  —  I  £  c^x  -f^l  =  — (c'‘x-f^).  (2.35) 

e  V  I  e'v 


Since  *  has  been  assumed  to  be  a  minimizer  of  the  linear  programming 
problem  (1.1)-(1.4)  and  the  objective  function  is  normalized,  we  have 


c’’x  =  0.  (2.36) 

Therefore  the  objective  function  assumes  a  negative  value  at  u.  and  this  is 
absurd. 

Let  us  assume  now  that  xS  A  and  that  E(]uation  (2.21)  holds.  We  will 
show  that  X  =  A’e  is  a  minimizer  of  the  linear  programming  problem  with 
normalized  objective  function  (l.l)-(1.4).  In  fact,  from  (2.21)  we  have 


2,-(A“''-c)  =0.  (2.37) 

Using  A  instead  of  A  as  in  (2.10),  when  AA"  '  ^  is  of  rank  less  than  m  we 
have 


0  =  e^A'-‘'^^n„,.:,.(A<"M 

=  e'"A '  -“/2[  /  -  X^  'Wi  .A.V"A'' )  ■  ’  AA"  "*]  A"  '-c 

=  t^Xc-e^XA^{AX'’A^)'\\X’’c  =  e^’Ac  =  cV  (2.38) 

Therefore  we  have  that  x  is  a  feasible  point  where  c^x  =  0,  that  is,  x  is  a 
toinimizer  for  the  linear  programming  problem  with  normalized  objective 
hinction  (1.1)-(1.4).  ■ 
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Let  2  be  the  set  given  by 

2-{xefl'Ux-0.e'^x-l}.  (2.39) 

Lemm.v  2.8.  Let  x*  be  a  nondegenerate  minimizer  of  tite  linear  program¬ 
ming  problem  with  normalized  objective  Junction  (1.1)-(1.4).  and  let  h(x)  be 
given  by  (2.11).  Then  we  have 

h(x*)=0.  (2  40) 

Moreover. 

Ah(x)  =  0.  xe2.  (2.41) 

e''h(x)=0.  xe2,  (2.42) 

where  2  is  given  by  (2.39). 

Proof.  In  fact  for  x  £  2  we  have  A.Ve  =  0.  e^.Ve  =  1.  and 
.AXrinYi.l.Yc)  =  0.  so  that 

.Ah(x)  =  -  .A.\n,,x, -(•Vc)  +  .A.\-ee’'.Vn, ,,,.(. Vc)  =0  (2.43) 

and 

e^h(x)  =  -e^X•^,,,,,-(Xc)+e^Xee^V^,,,v,-(-Xc)  =0.  (2.44) 

Moreoser.  from  Theorem  2.7  we  have  (2.40).  ■ 

Let  X  £  R".  and  £,  £/?'"“'  be  the  matrix  given  by 

£,=  DiaK(n,,„.(.Xc)).  x£fl"  (2  45) 

Let  /|,(x)£  fl'”'"  be  the  following  matrix: 

A{x)  =  -  ( /  -  .YeeMA  n,,,...V-'£,  -  [e^V H, Xc)]/.  x  £  R\ 

(2  46) 
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For  X  S  R"  in  (2.46)  we  will  use  A  instead  of  A.  as  in  (2.10).  when  .A.V  is  of 
rank  less  than  m.  .An  elementarv  computation  shows  that  the  matrix 
can  Ik*  defined  for  x  S  R“  so  that  is  defined  for  x  £  R". 
Let 


.y  =  {x  £  R"|.AA''A’^  is  invertihle} .  (2. 47) 

For  X  £  >,  Ji,(x)  is  the  jacohian  matrix  of  h(x)  with  respect  to  x.  Morc-iwer 
let  S(x*.p*)  be  the  ofien  ball  of  Lemma  2..5;  we  obsene  that  for  x  £  F  U 
S(x*.p*).  since  .AA'  is  of  rank  in.  the  matrix  is  well  defined 

o' 

and  cuntinticHis.  So  Jh(x)  is  continuous  for  x  £  F  U  S(x*.  p*),  and  since 
Fl,^ y., .(.V*c)  =  0.  we  have 


(2  4.S) 

From  Lemma  2.S.  we  conclude  that  any  solution  x*  of  the  linear  prutfrum- 
ming  problem  with  normalized  objective  function  (!.!)-( 1.4)  is  an  eijuilib- 
rium  point  of  (2.1.3).  that  is.  h(.x*)=  0. 

Hovveser.  due  to  the  sinsular  Jacobian  of  h(x)  at  x*  (that  is.  to  (2.4b)l.  the 
use  of  a  linearis  implicit  .A-stable  method  to  intecrate  the  initial  value 
problem  (2.13).  (2.14).  as  suggested  in  (8)  in  the  context  of  nonlinear 
programming,  will  not  produce  a  quadratic-ally  convergent  method  for  linear 
programming.  To  overcome  this  difficulty  we  introduce  a  new-  vector  field 
Rx)  SR”  defined  for  x  £  fl  "  given  by 

f(x)  = -(/- Aee’’)[xc-X.A^(.AXA'’)''.A.Vc].  x£fi".  (2.49) 

where  we  use  A  instead  of  A,  as  in  (2.10).  if  AXA^  .A  is  of  rank  less  than  m. 

Let  us  (xinsider  Rx)  for  x  £  F  U  S(x*.p*).  We  observe  that  .A-LA^  is 
invertible.  From  (2.49)  we  have  that  Rx)  is  a  continuousb  differentiable 
function  of  x  for  x  £  F  U  S(x*.p*).  For  later  purposes  we  observe  that  Rx)  for 
X  £  F  can  be  rewritten  as  follows; 

f(x)  =-(/- .\'ee'’)A-'^*n,.,.v. V'^-c),  x£F.  (2.50) 
or 

f(x)=  -A'/^(/-.V'-'-ee''X'''-)g(x.l).  x£F.  (2.51) 
From  Equations  (2.20).  (2.51)  and  Thc*orem  2.7  it  follows  that  if  x*  is  a 
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minimizer  of  the  linear  programming  problem,  then  R«*)  =  h(x*)^“  0,  that  is, 
X*  is  an  equilibrium  point  of  the  vector  fields  Wx).  Rx).  For  x  €  A  the  vector 
field  Rx)  can  be  obtained  as  the  steepest  descent  vector  associated  to  the 
function  c^x  with  respect  to  a  particular  metric.  In  [131  D.  A.  Bayer  and  J.  C. 
Lagarias  have  introduced  the  idea  of  looking  at  Karmarkar's  vector  field  h(x) 
in  teims  of  steepest  descent  directions. 

Let  Xg  be  a  feasible  point  of  the  linear  programming  problem  (1.1  Ml. 4), 
and  Fg  be  the  affine  subspace 

=  *u  ^  =0}.  (2.52) 

Lemm.4  2.9.  The  vector  field  Rx)  given  by  (2.51)  is  the  steepest  descent 
vector  associated  to  the  objective  function  b(x)  =  c^x  of  the  linear  program¬ 
ming  problem  (1.1)-(1.4)  restricted  to  Fg  O  F  with  respect  to  the  Riemanniun 
metric  G(x)=  =  Diag(x"*),  defined  on  the  positive  orihant  F.  where  F„ 

is  given  by  (2.52). 

O 

Proof.  We  consider  the  following  transformation  for  x  €  Fg  O  F: 

x  =  G-'-^(x)y  =  .\-'S.  (2.5.3) 

We  have 

/^(x(y))-(.Y'  -c‘'y.  (2 -...) 

and  Fg  assumes  the  following  form; 

Fg  =  Xg  -r  {u  €  =  0.  e^V'''-u  =  0)  (2.55) 

The  gradient  vector  of  f;(x(>  ))  with  respect  to  y  is 

db 

— -X'^-c.  (2.56) 


The  gradient  vector  db/dy  projected  on 


is  given  by 
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where  we  use  A  instead  of  A,  as  in  (2.10),  if  AXA^  is  of  rank  less  than  m. 
Since  .A.Ve  =  0  from  (2.57),  using  (2.7)  we  have 


I  =  A,K.4’^)‘‘.A.\'c  -  X'^-ee^Vc. 


(2.58) 


Since  e^XA’^(AXA^)" ‘AXc  =  o,  we  have 

5  =  (/-X'''-ee'-X’/-)n,.,^..,-(X''^c).  (2.59) 

Finally,  appKing  (2.53)  to  4,  we  have  that  the  gradient  vector  ^(x)  is  given 
by 


.(x)  =  .V''-4  =  .V'  -(/-X 


X'  -)n,,,,  -c),  (2.60) 


This  concludes  the  proof 
Let  \  be  given  bv  (1.8). 


■ 


Levivt.v  2.10.  Let  x  6  A,  and  r  be  the  rank  of  the  matrix  AXA^.  Then 
n r  . .. , .n  ‘  ( X ' /-c)  =  ( /  -  X ' ee^Y ' ^^ )  H, X '  ) ,  ( 2 .6 1 ) 


.A.V 


where  we  use  A  instead  of  A.  as  in  (2.10).  wher  r  is  less  than  m. 

ax'  *1 


Proof.  Let  0  <  r  <  m.  The  projector  Flj  ,  1  -  is  defined  as  follows: 


n 


•  =  /- 

■.AX'^^' 

■.4X‘/*  ■ 

’  .AX  ‘^-  ' 

T 

-AX'^^' 

AX^  ‘ 

Ic^X'  = 

.e'X'^V 

1 

erx</J 

.e^Y'^N 

) 

.  (2.62) 


Us  compute  the  matrix 
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Since  AXe  =  0  and  e^Xe  =  1,  we  have 


A/  = 

(.AXA’’)'' 

0 

0^  i 

1 

(2.63) 


An  elemental^'  computation  gives  us 

Rf  1  -  (  X  '/-c)  =  X  '^-c  -  .A'lA'’)  ‘  '.AVc  -  X  '"-ee^Vc.  (2.6A) 


From  e^XA^  =  0  we  have 


q  =  X'^-ee'’.\CA’’(  AXA^)''.A\’c  =  0  (2.&5) 

and 


n 


AV 


. ,  -(X'/-’c)  =  n,,,.  -c)  -  X'^^ee’-X 


»' V' 


1/2 


c  +  q,  (2.66) 


With  an  eas>  computation  from  (2.66)  we  obtain  (2.61).  Let  r  =  0.  and 
O  e  fi " ' "  be  the  null  matri.K.  We  hav  e  that 


O 

er_^./2 


=  (e*^.\''  ■)'  and  n„-=/. 


so  that  (2.61)  holds. 


LeM\iA  2.11.  Let  x*  be  a'nunde^cnerate  miniinizer  of  the  linear  pro¬ 
gramming  problem  uith  normalized  objectiic  function  (l.l)-(l  4);  let  fix)  be 
given  by  (2.49)  and  1  be  given  by  (2.39).  Then  we  have 

f(x*)=0  (2.67) 

moreover 

.Af(x)  =  0.  xel.  (2.68) 

e’’f(x)  =  0.  xel.  (2.69) 


(2.79) 
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Proof.  From  Lemma  2.11  we  have 

dx 

A— -Af(i)=0,  (2.80) 

dx 

e^— =  e^f(*)  =0,  (2.81) 

so  that  the  thesis  follows  immediately  from  the  .issumption  (2.74),  (2.75)  on 

Xg  and  the  fundamental  theorem  of  calculus.  B 

For  *  e  fi*  let  £  G  fl"*"  be  the  matrix  given  by 

£  =  Diig{A^AXA^)''AXc),  x  6  R\  (2.82) 

and  C  e  £"*"  be  the  matrix  C  =  Diag(c).  Let  J(x)e  fi"*"  be  the  following 
matrix; 

J(x)  =  -(/-.’L4''(.A.X.A^)‘'A-Xee’')(C-  £)+{e’'.Vc)/,  x  e  fl". 

(2.a3) 

where  we  use  -A  instead  of  .A,  as  in  (2.10),  if  AXA^  is  of  rank  less  than  m. 
Let  .y'  be  the  set  (2.47).  .An  elementary  computation  shows  that  for  x  e  /F", 
}(x)  is  the  Jacobian  matrix  of  Rx)  with  respect  to  x.  Moreover  let  S(x*,p*)  be 
the  open  ball  of  Lemma  2.5.  We  observe  that  for  x  G  F  U  S(x*.p*).  since  the 
matrix  AXA^  is  invertible,  /(x)  is  continuous  for  x  G  P  U  S(x*.p*). 


Theorem  2.13.  Let  us  assume  that  the  linear  programming  problem 
(1.1)-(1.4)  has  a  unique  nondegenerate  minimizer  x*,  and  let  Jlx*)  be  given 
ly  (2.83).  Then  ]{x*)  is  invertible  as  an  operator  restricted  to  the  subspace 

—1  .  That  is.  /(x*)v  *  0  for  each  v  #  0  such  that 
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Proof.  First  of  all  we  show  that  for  i  =  1.2 . n  we  have 

( C  -  £)„  =  0  if  and  only  if  i  *  #  0. 

Let  X*  =■  Diai»(x*).  From  Theorem  2.7  fiir  a  =  1  we  have 

X*'^-(C-£)  =  Diai;((n,,v.-,-(.V‘'-c)))  =  0.  (2.S4) 

Since  C  —  E  is  a  diagonal  matri.x.  from  (2.S4)  we  have  that  X*  *  0  implies 

(C—  £)„=0  for  i  =  1.2 . n.  Let  us  show  that  (C— £)„=0  implies 

X*  *  0  for  i  =  1,2 . n.  In  fact  if  we  assume  that  there  e.vists  h  such  that 

(C  —  £);,/,  =  0  and  =  0.  then  from  the  assumption  that  x*  is  a  nondegen- 
erale  minimizer  of  the  linear  programming  problem  (l.l)-(l.-l)  it  follows  that 
there  exists  a  pivot  transformation  that  makes  the  hth  component  of  ** 
nonzero.  Let  y*  be  this  new  basic  feasible  solution  corresponding  to  x*  via 
the  pivot  transformation.  Since  only  one  pivot  operation  has  been  made,  the 
nonbasic  components  other  than  the  /ith  component  are  still  nonbasic.  that 

is,  (/•  »0  for  each  such  that  i*-0.  Let  )*=  Diag(y*.y* . tf*). 

From  (2.84)  we  have 

y*(C-£)=0  (2.85) 

Moreover,  since  e^Y*A^  =  0,  from  (2.85)  we  have 

0  =  e^Y*{C-  £)e  =  e^Y*Ce-e^Y\V{AX\V)'' .\X*c  =  e’^Vc  =  c'^yV 

(2,86) 

Therefore  y*  would  be  a  new  mininiizer  of  the  linear  programming  problem 
(11)-(1.4),  different  from  x*.  and  this  is  absurd. 

We  have  /(x*)*  =(C— £)*.  In  fact  it  is  obvious  that  vE{C—E)^ 
implies  v  e  y(x*)*.  Moreover,  let 

.M  =  [A'’(.AX*A'')''.A-ee"](C-  £).  (2.87) 

Then 

;(x*)-(C-£)-X*.M.  (2.88) 

so  that  /(x*K  =  0  implies  (C  —  £)v  =  X*M\  Since  X*  is  a  diagonal  matrix. 
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we  obtain  {X*Mv),  =  (X*)„(.Wv),.  i  =  1.2 . n.  VVe  have  two  cases; 

(i)  X*  =  0,  which  implies  ((C  -  EM,  =  0: 

(ii)  X*  #  0,  which  implies  (C  -  E)„  »  0. 

Summarizing,  we  have  that  ((C-EM, —  0,  i»»l,2 . n.  that  is.  v  e 

(C-E)^. 

Now  let  u  €  fl"  be  such  that 


Au  =  0.  e’’u  =  0. 


(2.89) 


We  assume  that  u  e  ](%*)  ~ ;  since  x*  e  (C  —  E)  *  ,  then  z  =  i*  +  u  e  J(x*)^ . 
Moreover, 


Az  =  0.  e’^z  =  1.  (2.90) 

and  z  e  implies  z  6(C  -  E)^.  If  z  e  (C  -  E)"^,  then  c,  =  0  for  each 

i  such  that  i*  =  0;  this  condition,  together  with  (2.90),  is  a  characterization 
of  the  minimizer  x*  of  the  linear  programming  problem  (1,1)-(1.4).  There¬ 
fore  u  =  0.  This  concludes  the  proof.  ■ 


Theorem  2.14.  Let  x*  be  the  unique  nundeitcnerate  minimizer  of  the 
linear  programming  problem  with  normalized  objective  function  (1.1)-(1.4). 
and  Rx)  be  given  by  (2.49).  Hr  consider  the  initial  value  problem 


dx 

—  =  f(x).  (2.91) 


1 

.x(0)  =  -e.  (2.92) 

n 


Then  a  solution  x(f,(l/n)e)  of  (2.91).  (2.92)  e.Tisfs  for  t  e(0.  =c),  and 

lim  x|t.— ej  =  X*  (29.3) 


Proof.  The  standard  existence  and  unirjueness  theorems  for  the  initial 
value  problem  for  ordinars  differential  equations  guarantee  that  the  solution 
of  (2.91).  (2.92)  exists  locallv.  From  Lemma  2.9  it  follows  that  Rx)  is 
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tangential  to  5A,  so  that  from  Lemma  2.12  and  the  fact  that  (l/;i)e  €  A  we 
have  that  x(<,(l/n)e)e  A.  Moreover  for  *  €  we  have 


d 

7t 


(2.5«) 


that  is.  the  objective  function  c^x  is  monotonically  decreasing  along  the 
trajectorv  xit.(l/n)e).  Since  the  minimun)  of  on  is  zero,  x*  is  the 
unique  minimizer  of  c^x  on  .\.  and  fix*)  =  0,  from  C.  Sansone  and  R.  Conti 
l21.  p.  .31]  we  have  that  x*  is  the  uni(]ue  limit  point  of  x</.(l/n)e)  and 


(29S) 


This  c-oncludes  the  proof 


3.  THE  QUADRATIC  ALGORITHM  FOR  LINEAR  PROGRAMMING 


_^t  xe  R".  D  C.  R"  be  an  open  set.  and  D  be  the  closure  of  D:  let 
w:DcR"-«R"  be  a  function  continuously  differentiable  in  D.  whose 
Jacobian  matrix  is  denoted  by  Q(x)  =  dw/dx.  Let  us  consider  the  initial 
value  problem 


—  =  w(x).  (3.1) 

x(0)  =  x„.  *„eD.  (3.2) 


Let  /  be  the  n  X  n  identity  matrix,  h^  >  0.  k  =  0.1.2 .  be  a  sequence  of 

Mepsizes.  and  Then  any  solution  x(tj)  of  (3.1).  (3.2)  can  be 

approximated  with  x*  computed  as  follows; 


(/-/.tQ(x‘)ls* 


*0' 

(33) 

/I,w(x‘).  0.1.2 . 

(3.4) 

x*+$*.  it -0.1. 2 . 

(3.5) 
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The  numerical  scheme  (3.3)-(3.5)  to  integrate  the  initial  value  problem  (3.1), 
(3.2)  is  A-stable  and  linearly  implicit,  and  has  Jieen  studied  by  J.  D.  Lambert 
and  S.  T.  Sigurdsson  in  [22]. 

Let  Rx)  be  the  vector  field  given  by  (2.49),  and  J(x)  be  its  Jacobian 
matrix  given  by  (2.83).  We  will  apply  the  numerical  scheme  (3.3)-(3.5)  to 
the  initial  value  problem 

dx 

-^  =  f(t).  (3.6) 

1 

x(0)  =  -e  (3.7) 

n 

considered  in  Section  2.  Let  V  be  given  by  (2  9).  By  Lemma  2.5  there  e.xists 
a  neighborhood  S(x*.p*)  of  x*  such  that  fix)  is  a  continuously  differentiable 
function  in  f  uS(x*,p*). 

LtMMA  3  1.  Let  X*  be  the  unique  minimizer  of  the  linear  programming 
problem  with  normalized  objective  function  (1.1)-(1.4).  Let  us  apply  the 
numerical  scheme  (3.3)-(3..5)  to  the  initial  value  problem  (3.6).  (3.7).  More¬ 
over  let 


/ij  *  (c%c)  for 

L  =0.1.2 . 

(3.S) 

1  —  h^J{x'‘)  be  invertible 

for  A- =0.1.2 . 

(3  9) 

Then  the  sequence  {x*  *  ').  k  =  0. 1.2 .  generated  by  (3.3)-(3.5)  exists,  and 

s*  satisfies 

As‘  =  0.  A=0.1.2 .  (3.10) 

eV  =  0.  A  =0.1. 2 .  (3.11) 

Proof.  We  note  that  x'’  =  (l/n)e  is  a  feasible  point  of  the  linear 
programming  problem  (1.1)-(1.4)  and  that  is  defined  by 

[l-hj(x‘‘)]s^~hj{x'‘).  11  =  0.1.2 .  (3  12) 

In  Lemma  2.11  it  has  been  showu  that  if  x*  is  a  feasible  point  for  the  linear 
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programming  problem  {1.1)-(1.4),  we  have 

Af(x*)=0. 

=0. 

so  that  applying  A  to  both  sides  of  (3.12).  we  obtain 

A[l-hj{x^)]s^  =  0. 

From  (2.83)  we  have 


Therefore  we  have 


A;(x‘)  =  (e%c)A. 


l-/it(e%c)lA$‘  =0. 


so  that  from  (3.8)  we  have  As*  =  0,  k  =  0, 1,2 . 

Moreover  it  is  eas>’  to  verifr  that 

e";(x*)  =  (e%c)e"/. 

so  that  from  Lemma  2. 1 1  we  have 

e"(/-A»/(x*)]s‘=0. 

"'hich  implies 


[l-/.4(e%c)]eV=0. 

so  that  from  (3.8)  we  have  e^s*  =  0.  A:  =  0. 1.2 . 

Let  DqR"  be  an  open  set  and  Dq  C  D  be  a  convex  set. 


(3.13) 

(3.14) 

(3.15) 

(316) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 


DcFiNtTio.N  3.2.  Let  w(x).  DS  R"  be  a  continuously  differen- 

function.  Let  {  e  fl",  be  an  open  neighborhood  of  4.  and  T ;  D  X 
C  fi"  X  /I”  -»  L(fl*),  where  L(fl")  is  the  set  of  the  n  X  n  matrices.  Then 
^^*-4)  is  a  consistent  approximation  to  the  Jacobian  matrix  Qif)  of  w(x)  on 
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Dq  C  O  if  0  €  if"  is  a  limit  point  of  D^  and 

^lim  r(x.|)=P(x)  (3il) 

uniformly  for  x  e  Dq.  Moreover,  if  there  exist  two  constants  c  >  0  and  r  >  0 
such  that 

l|Q(x)-r(x.4)|l<c||4ll  (3.22) 

for  each  x  e  D„  and  5^D^nS(0.  r).  where  S(0.  r)  =  {|  e  fl’"  |  II^H  <  r). 
then  T(x.  4)  is  a  strongly  consistent  approximation  to  Q(x)  on  D„. 

Lemma  3.3.  Let  DcR"  be  an  open  set.  and  w;Dcfl"-»i?"  be  a 
continuottsly  differentiable  function  on  the  convex  set  D„  C  D.  Let  Q(x)e 
Lip^(  D|,).  that  is.  let  Qix)  be  a  Upschitz  continuous  function  for  x  €  Dy  with 
Upschitz  constant  y  >  0.  Then 


II  w(y)  -w(x)-  P(x)(y-x)||<  ^llx-yil’ 


(3.23) 


for  each  x.  y  S  D„. 


Proof. 
p.  73]. 


See  J  M  Ortega  and  W.  C.  Rheinboldt  [20.  Theorem  3.2.12. 


TiitoBEM  3.4.  Let  X*  be  the  unique  nondegenerate  ininiinizer  of  the 
linear  programming  problem  with  normalized  objective  function  (1  1)-(1.4). 
and  let 


8  = 


.A 


and  B  ~  be  given  by  (2  S).  Moreover,  let  fix)  be  given  by  (2.49).  and 


K 


it  -0.1.2 . 


(3.24) 


where  >  a  >  0  is  a  bounded  sequence 
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there  exists  p|  >  0  and  an  open  neighhorhuod  S(x*.p|)  =  {x  e  fl"l|lx-x*||  < 
p,)  of  X*  such  that  if  4.1  ^  S<x*.p,)n  S.  uhere  S  is  given  by  (2.39).  then  the 
sequence  {x*},  it  =  0. 1,2 .  generated  by 


=  (3.2.5) 

[/-/iiy(xM](x‘*'-x‘)  =  /ijf(x‘).  it=0.1.2 .  (326) 

where  }(x)  is  given  by  (2.83)  and  the  linear  system  (.3.26)  is  solved  in  B  ~ .  is 

well  defined.  x**'eS(x*.p,)nl  for  k  =  0. 1.2 .  and  {x^}.  k  =  0.1.2 . 

converges  quadratically  to  x*. 


Proof.  We  armif  b\  induction  on  the  index  k.  Let 

1 

^^  =  -.  i=0.1.2 .  (3.27) 


and 


<f(x‘.fi)= -Ii/  +  y(x‘).  it  =0,1. 2 .  (3.28) 

We  obsene  that  (3.26)  can  be  rexxritten  as  follows; 

-<t>(x‘.^J(x*'' -x‘)  =  f(x‘),  it  =0.1. 2 .  (3.29) 

It  is  easy  to  see  that  <t>(x,^)  is  a  strongly  consistent  approximation  of  Jlx)  on 
r  when  f  goes  to  zero.  We  have  seen  in  Lemma  3. 1  that  if  v  6  fl "  and 
**  S  i  we  have 


A<t>(x‘.|i)v  =  (e%c-44)Av  (3.30) 

e’’4>(*‘.4^)v  =  (e^X■^c-^^)e'■v.  (3.31) 

From  (3.24)  and  (3.27)  it  follows  that  4*  *  e^X^c  for  t  =  0. 1,2 .  so  that 

have 


‘J>(x*.ft)v  €  B  ■  if  and  onl\  if  v  e  B  . 


(3.32) 
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From  Theorem  2.13  we  have  that  J(.x*)  restricted  to  B  is  invertible  and 

||(y(x*)li;^'(/(**)  -  oil,.  I  <  1  (3-33) 

for  x*£S(x*,p,)  in  a  suitable  neighborhood  of  x*  and  in  a  suitable 
neighborhood  U  of  zero.  The  perturbation  lemma  (see  J.  M.  Ortega  and 
W.  C.  Rheinboldt  [20.  Lemma  2.3.2.  p.  -IS])  implies  that  the  inverse  of  the 
linear  operator  <IKx*.^4)  restricted  to  the  subspace  B"^  e.tists  when  x*  £ 
S(x*,p,),  £  t/.  From  Lemma  2.11  we  have  that  x*  £  S  implies  Rx*)  £  B 

When  X*  £  S(x*.p,)n  1  and  from  the  fact  that  Rx‘)£B^  and 

(3.32),  (3.33)  it  follows  that  x**'  is  well  defined  and  x*'''£l.  Moreover 
there  exists  ij  >  0  such  that 

||‘l’(*.^)l8*.'||<  n-  x£S(x*,p,),  f£l’.  (3.34) 

and  we  have 

IU“'-x*ll  =  ||<l>(x‘,^*)|;.'[4>(**.fJ(x‘-%‘)-f(x*)]j 

<’j(ll*(*‘.^*)-y(*‘)||+||;(x*)-/(x‘)!i)!ix‘-.x*ii 

+  »,l|f(x*)-f(x‘)-y{x*)(x‘-x*)||.  (3.35) 

Since  we  can  always  choose  p,  >  0  such  that  y(x)£  Lip^(S(x*.p,))  for  some 
y  >  0.  from  Lemma  3.3  we  have 

ny 

llx‘’  '  -x*ll  <  f*r,|lx‘  -  x*ll+  7,yllx‘  -  x*!!*  +  y  llx*  -x*|l' 

<w(x‘.fj!!x‘ -.x*ll.  (3.36) 

where 

w(x‘.fi)  =  +  ;»7yllx* -x*||.  (3.37) 

The  neighborhoods  S(x*.p,)  and  C  can  be  chosen  in  such  a  way  that 

w(x*.fj<x<  1.  X*  £S(x‘.p,)nS.  (3.38) 
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From  (3.36)  and  (3.38)  we  have 

llx‘*' -x‘IUxllx‘-x*l|.  x*eS(xVp,)nS.  (3.39) 

Therefore  S  S(x*.P|)0  1.  Moreover  we  ha\e 

-x‘li<x*-'llx"-x*||.  (3.40) 

so  that  the  iterates  (x*),  k  =0,1.2 .  given  by  (3.26)  are  well  defined  if 

x"  6  S(x*,P|  )n  1  and 

limx‘=x*.  (3  41) 

Moreover,  since  Rx)  is  a  continiiousK  differentiable  fiinction  in  a  neighbor¬ 
hood  of  X*  and  Rx*)  =  0.  there  e.xists  a  constant  M  >  0  such  that 

l|f(x‘)||<.'/!lx‘ -X*||.  X*  eS(.v*.p,)nS.  (3,42) 

From  (3.24)  we  obtain 

Ilx‘-X*ll,  x*eS(x*.p,)n^.  (3.43) 

w 

That  is.  the  sequence  (x*)  converges  quadraticalls  to  x*.  This  concludes  the 
proof.  ■ 


4.  NLMERIC.\L  EXPERIMENTS 


V\'e  begin  by  comparing  the  computaMonal  cost  of  a  step  of  the  algorithm 
introduced  in  Section  3  with  that  of  a  step  of  the  simplex  algorithm  or  of  a 
^iep  of  Karmarkar's  algorithm. 

We  consider  the  linear  programming  problem  in  the  canonical  form 
(1.1)-(1.4).  It  is  easy  to  veriR  that  the  computational  cost  of  a  step  of  the 
simplex  algorithm  is  given  by 


mn  +  lower  order  terms . 


(4.1) 
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The  computationiil  cost  of  one  step  of  Karmarkar  s  algorithm  is  essentially 
due  to  the  computation  of  the  matrix 

AX^A^  (4.2) 

and  to  the  solution  of  the  m  X  m  linear  system 


( AXW)y  =  .AX-c. 


(4.3) 


Since  the  matrix  AX\\^  is  symmetric,  its  (imputation  re<juires 


m'n 

^  -  +  lower  order  terms  (4.4) 


elementary  operations,  while  the  solution  of  the  linear  system  (i  Z)  requires 


m' 


—  +  lower  order  terms 

6 


(4.3) 


elementary  operations.  Since  m  <  n.  we  can  conclude  that  the  computational 
cost  of  one  step  of  Karmarkar  s  algorithm  is  roughly 


jn'’ +  lower  order  terms  (4  .6) 

elementary  operations.  This  cost  can  be  reduceJ  using  some  special  proce¬ 
dures;  for  example,  in  [2]  Karmarkar  has  shown  that  the  use  of  successise 
rank-one  modifications  to  compute  f4.2).  (4.3)  reduces  the  "asera^e’  compu¬ 
tational  cost  of  each  step  to 

cn~^  +  lower  order  terms  (4.7) 


elementary  operations,  for  some  constant  c  >  0. 

The  computational  cost  of  one  step  of  the  algorithm  introduced  in  Section 
3  is  essentially  due  to  the  computation  of  the  matrix 


.VUJ. 


(4.3) 
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to  the  solution  of  the  linear  sssteni 


(  =  .-EVc.  (-J.9) 

to  the  eoniputation  of  the  matrix 

.V(AXA- ,  'a  (4,10) 

that  appears  in  the  expression  of  the  Jacobian  J(x)  (2. S3),  and  to  the  solution 
of  the  linear  sxstem  (3.26).  The  coniputational  ct)sts  of  (4,S).  (4.9)  are 
analogous  to  those  of  (4.2).  (4.3)  respectixeK.  .Moreover  the  computational 
cost  of  the  solution  of  the  linear  system  (3.26)  is 


n  ’ 

—  -i-  lower  order  tenns 


(4.11) 


elementary  operations. 

In  order  to  compute  (4.10)  we  use  the  Chi  lesky  decomposition  of 
that  is. 


AX.\^  =  UJ. 


(4.12) 


where  L  €  fl'"’'"*  a  nonsingular  lower  triangular  matrix.  So  we  t  ive 


(.A.\LA^)’'=(L'')"l- 


(4  13) 


anc. 


a^(ax.^^)''a  =  {l-'a)^{L-'a). 


(4.14) 


-ince  in  order  to  compute  L‘\\ 


m  ’Ti’n 

—  + - +  lower  order  terms 

6  2 


(4.15) 


^It.nenUrv  operations  are  neeessarv.  and  the  matrix  (L"U)^(L“‘A)  is 
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symmetric,  the  compMitatiunal  cost  of  (4.10)  is 


m^n 


(4.16) 


elementary  operations.  Since  m  <  n,  we  can  conclude  that  the  computational 
cost  of  one  step  of  the  algorithm  introduced  in  Section  .3  is 

+  lower  order  terms  (4.17) 


elementary  operations.  .Moreover,  if  we  use  successive  rank-one  mixlifica- 
tions,  as  proposed  in  [2],  we  can  decrease  the  "average"  computational  cost  of 
each  step  to 


c'n*  *  + lower  order  terms  (4. IS) 

elementan  operations,  for  some  constant  c'>0.  Moreover,  to  improve  the 
value  of  c'  it  is  possible  to  use  any  combination  of  the  ideas  proposed  in  [7. 
23,  24]. 

To  conclude,  the  computational  cost  of  one  step  of  the  algorithm  intro¬ 
duced  in  Section  3  is  of  the  same  order  as  that  of  one  step  of  Karmarkar  s 
algorithm,  while  one  step  of  the  simplex  algorithm  is  much  cheaper.  How¬ 
ever.  due  to  the  quadratic  convergence  of  our  algorithm,  we  expect  that  the 
number  of  iterations  needed  to  solve  a  linear  programming  problem  to  a 
given  accuracy  should  be  approximately  independent  of  the  problem  size  n. 

We  present  now  some  numerical  results  that  support  our  expectation. 
The  algorithm  described  in  Section  3  has  been  implemented  using  two 
sptecial  expedients  to  avoid  failure  due  to  the  ill-conditioning  of  the  problem 
considered. 

The  matrix  .A  given  by  (1.2)  is  replaced  with  the  matrix  .AS 

jq  reduce  its  condition  number;  .A  is  obtained  usine  the  singular  value 
decomposition  of  .A.A^  This  decomposition  has  a  computational  ctist  of  order 
n^.  so  it  costs  the  same  as  one  step  of  the  algorithm  described  in  Section  3 
Let  =  Q^DV  be  the  singular  value  decomposition  of  .\A^\  then  the 
matrix  .A  is  given  bv 


A  =  {Q^D*V).\.  (4  19) 

where  D*  G  a  diagonal  matrix  such  that  ( D*)„  =  1/D„  if  D„  >  0 

and  ( £)*)„  -  1  if  D„  “  0  for  i  =•  1 . m. 
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In  the  first  steps  of  our  algorithm  (i,  <  5  in  our  numerical  experi¬ 
ments),  the  Riemannian  metric 

C(x)  =  A-'  =  Diag(x-')  (4.20) 

is  replaced  with 

C(x)  =  M,X-\  (4.21) 

where  =  Diaglx^ ')  and  Xj  is  the  current  point  at  step  k.  We  note 

that  in  order  to  apply  our  algorithm  is  not  necessary  to  have  a  normalized 
objective  function — that  is.  is  m)t  necessar\-  to  know  the  value  of  the 
objective  function  at  the  minimizer.  However,  in  our  numerical  experiments 
we  use  test  problems  with  normalized  objective  function.  The  stopping  rule 
used  is 


=c''x‘ «  1  0xl0-'’jc’’^j.  (4  22) 

We  have  considered  ten  test  problems.  Problem  1  ( Z  I R 1 )  is  a  problem 
coming  from  the  operation  of  an  industrial  plant  in  central  Italy.  The  other 
problems  come  from  the  System  Optimization  Laborators  at  Stanford  Uni¬ 
versity  and  have  been  made  available  to  us  through  netlib  (25).  The 
numbers  of  variables  (n)  and  of  constraints  (m),  shown  in  Table  1,  are  those 
relative  to  the  problems  in  canonical  form.  Finally,  it  denotes  the  index  of  the 
first  step  that  verifies  the  stopping  rule  (4  22). 

We  note  that  in  Table  I.  while  n.m  varv-  by  an  order  of  magnitude,  the 
number  k  of  steps  needed  to  solve  the  problem  varies  onlv’  by  a  factor  of  two. 
Moreover,  test  problems  with  n.m  are  solved  in  about  ten  steps. 


TABLE  1 


Test  problem 

m 

n 

k 

*'4 

1. 

ZIR1 

304 

543 

21 

2.41D-10 

2. 

XOLITTLE 

57 

141 

21 

3.16O-09 

3. 

AFIRO 

2S 

54 

12 

1  52ol2 

4. 

BEACONFD 

173 

298 

20 

3.91  d-09 

5 

BLEND 

75 

117 

21 

1.470-12 

6. 

ISRAEL 

175 

319 

17 

1940-10 

t  . 

SC10S 

106 

166 

13 

1  36i>ll 

8. 

SC50A 

51 

81 

14 

1.480-14 

9. 

SCSOB 

51 

81 

11 

7.84  o- 10 

10. 

SHAREZB 

97 

167 

21 

1.780- 10 
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The  algorithm  has  been  cxxled  in  fortran  and  tested  on  a  VAX/VMS 
Version  V5.1  in  double  precision  arithmetic. 
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